TreatmentSimulator_SNandNC_Fire_ProductionNotebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c5e31164-d11e-45d6-98b4-b641623b8ddc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import os\n",
    "import math\n",
    "from datetime import datetime\n",
    "from src.utils import ee_treatments\n",
    "from src.utils import runs\n",
    "from src.utils.check_exists import check_exists\n",
    "from src.utils.yml_params import get_export_params\n",
    "yml_file = os.path.join(os.getcwd(),'config.yml')\n",
    "CRS,EXPORT_SCALE = get_export_params(yml_file)\n",
    "\n",
    "ee.Initialize(project='mas-gee')\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# In development phase\n",
    "tag = 'ee_treatments_production'\n",
    "ee.data.setWorkloadTag(tag)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2bc904bc-7607-4e66-a384-2fedde5094a8",
   "metadata": {},
   "source": [
    "## Million Acres Strategy - Treatment Simulator - **Testing Notebook**\n",
    "#### Creates Randomized Treatment Scenario Landscapes - exports as ee.Images\n",
    "#### This Notebook is for Run IDs: \n",
    "|RunID | CA_REGION | INTENSITY_ID | PRIORITY |\n",
    "|------|-----------|--------------|----------|\n",
    "| 3  | Sierra Nevada | 500k | Fire |\n",
    "| 6  | Sierra Nevada | 1m   | Fire |\n",
    "| 9  | Sierra Nevada | 2m   | Fire |\n",
    "| 12 | North Coast   | 500k | Fire |\n",
    "| 15 | North Coast   | 1m   | Fire |\n",
    "| 18 | North Coast   | 2m   | Fire |\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ba63a821",
   "metadata": {},
   "source": [
    "Run Settings: Change `RUNID` below to setup the analysis and output file names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e8b92d55",
   "metadata": {},
   "outputs": [],
   "source": [
    "RUNID = 'RunID18' "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b8cd85d3",
   "metadata": {},
   "source": [
    "Auto-calculate the necessary information for the Scenario\n",
    "\n",
    "Double-Check the Output Calcs are Correct against the [Spreadsheet](https://docs.google.com/spreadsheets/d/1Gnl6SO5kOkj4Ne1JdXzW4bp03824Zb1Cp_2HVuCY3Lk/edit#gid=0)\n",
    "\n",
    "(there may be small rounding errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b8dd49b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output EE Image Path: projects/mas-gee/assets/treatment_scenarios/RunID18/NC_Fire_2m_20230707\n",
      "TREATMENT SIZE (m^2): 404685.64224\n",
      "TREATMENT SIZE (Ac): 100\n",
      "Area to Treat Per Year (m^2): 2266239596.544\n",
      "Area to Treat Per Year (Ac): 560000\n",
      "Area to treat per 5-year interval (m^2): 11331197983\n",
      "Area to treat per 5-year interval (Ac): 2799996\n"
     ]
    }
   ],
   "source": [
    "# string identifiers for output file\n",
    "CA_REGION = runs.key[RUNID]['RegionID'] # NC = North Coast; SN = Sierra Nevada\n",
    "INTENSITY_ID = runs.key[RUNID]['IntensityID'] \n",
    "PRIORITY = runs.key[RUNID]['PriorityID'] \n",
    "\n",
    "# construct output path\n",
    "today_string = datetime.utcnow().strftime(\"%Y-%m-%d\").replace(\"-\", \"\")\n",
    "output_image = f\"projects/mas-gee/assets/treatment_scenarios/{RUNID}/{CA_REGION}_{PRIORITY}_{INTENSITY_ID}_{today_string}\"\n",
    "print(f'Output EE Image Path: {output_image}')\n",
    "\n",
    "# treatment size is static and set in runs.py\n",
    "treatment_size_ac = runs.treatment_size_ac \n",
    "TREATMENT_SIZE = runs.treatment_size_sqm \n",
    "print('TREATMENT SIZE (m^2):',TREATMENT_SIZE)\n",
    "print('TREATMENT SIZE (Ac):',treatment_size_ac)\n",
    "RADIUS = math.sqrt(TREATMENT_SIZE)/2  # in meters, rough square radius is A = (side/2) solve for side (side = sqrt(A) )\n",
    "\n",
    "sqm_to_ac_x = 0.000247105\n",
    "\n",
    "# pull prescribed treatment area in both units, we use sqm for analysis and acreage for reporting\n",
    "total_treated_area = runs.key[RUNID]['sqmPerYear']\n",
    "print(f\"Area to Treat Per Year (m^2): {total_treated_area}\")\n",
    "total_treated_area_ac = runs.key[RUNID]['acreagePerYear']\n",
    "print(f\"Area to Treat Per Year (Ac): {total_treated_area_ac}\")\n",
    "\n",
    "# we run the treatment generator as 5-year interval outputs (4 run intervals covering the 20-year period)\n",
    "# so we actually need to total_treated_area*5 to the ee_treatments() generator\n",
    "area_per_iteration = int(round(total_treated_area*5))\n",
    "print(f'Area to treat per 5-year interval (m^2): {area_per_iteration}')\n",
    "area_per_iteration_ac = int(round(area_per_iteration*sqm_to_ac_x))\n",
    "print(f'Area to treat per 5-year interval (Ac): {area_per_iteration_ac}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d06f1c00-707f-4c52-a9ea-1c9df7a9e419",
   "metadata": {},
   "source": [
    "Input Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6a94304a-b6c8-4a6b-b1b5-a2e4e504ae7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sierra Nevada(SN), South Coast(SC), Central Coast(CC), North Coast(NC)\n",
    "all_hucs = ee.FeatureCollection(\"projects/mas-gee/assets/TxPrctRankRrkWipRffc\")\n",
    "# Treatable Vegetation Constraint Layer\n",
    "# For reference Herb (1), Shrub (2), Hardwood (3), softwood (4), non treatable veg is masked\n",
    "veg_constraint = ee.Image(\"projects/mas-gee/assets/HerbShrubHardSoftImage\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "88dbf1b5",
   "metadata": {},
   "source": [
    "Specific Set up for North Coast and Sierra Nevada Fire Prioritization. \n",
    "\n",
    "see [Treatment Allocation Runs sheet](https://docs.google.com/spreadsheets/d/1Gnl6SO5kOkj4Ne1JdXzW4bp03824Zb1Cp_2HVuCY3Lk/edit#gid=0)\n",
    "\n",
    "* 1-25th percentile HUC group in Year Interval 1\n",
    "* 25th-50th percentile HUC group in Year Interval 2\n",
    "* 50-75th percentile HUC group in Year Interval 3\n",
    "* 75th-100th percentile HUC group in Year Interval 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8318bc6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "year_interval_ids = ee.List(['Y1to5',\n",
    "                             'Y6to10',\n",
    "                             'Y11to15',\n",
    "                             'Y16-20'\n",
    "                             ]\n",
    "                             )\n",
    "huc_filters = [ee.Filter.eq('TxBpPrc', 25),\n",
    "               ee.Filter.eq('TxBpPrc', 50),\n",
    "               ee.Filter.eq('TxBpPrc', 75),\n",
    "               ee.Filter.eq('TxBpPrc', 100)\n",
    "               ]\n",
    "# zipping ee.Lists is a good way to parallelize complicated things in GEE with .map()\n",
    "zipped = year_interval_ids.zip(huc_filters)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6fb01d8a",
   "metadata": {},
   "source": [
    "Pay attention to `ptOvershoot` value passed in `ee_treatments()` call. It is your dial. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fabc7615",
   "metadata": {},
   "outputs": [],
   "source": [
    "PT_OVERSHOOT = 2.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1956e6fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for QA test runs, get biggest and largest HUC poly per percentile group, may need to coded into this treat_by_year_huc() \n",
    "def treat_by_year_huc(l):\n",
    "    \"\"\"\n",
    "    this is to be map()ed across zipped (year intervals and percentile huc filters) to parallelize use of ee_treatments()\n",
    "        result is one ee.Image of treatment areas per yearInterval/HUC Group element in zipped     \n",
    "    \"\"\"\n",
    "    zipped_list = ee.List(l)\n",
    "    ranks_filter = ee.Filter(zipped_list.get(1))\n",
    "    \n",
    "    # Get HUC group\n",
    "    huc_group = (all_hucs.filter(ee.Filter.eq('RRK_Rgn', runs.key[RUNID]['Region Filter Name']))\n",
    "                .filter(ranks_filter)\n",
    "                # .limit(2) # for testing, reduce compute for development\n",
    "                ) \n",
    "    \n",
    "    # area_per_iteration is the prescribed area \n",
    "    # pixel_value_iteration is the yearInterval code to assign the treated pixels\n",
    "    pixel_value_iteration = ee.Number(ee.List(zipped).indexOf(l)).add(1)\n",
    "    \n",
    "    # this is for reshuffling the random treatment placements, may be useful for re-treatments of the same huc groups\n",
    "    # we would pass a unique number to 'seed' arg in ee_treatments() each time \n",
    "    # if you don't need to reshuffle the random seed each run, don't define 'seed' arg and function will use its default seed every time\n",
    "    seeds = ee.List([10,20,30,40])\n",
    "    seed = ee.Number(ee.List(seeds).get(ee.List(zipped).indexOf(l)))\n",
    "    \n",
    "    treated_area_img, properties = ee_treatments.ee_treatments(hucs=huc_group,\n",
    "                                                  prescription=area_per_iteration, # area_per_iteration_testing , reduce compute for testing\n",
    "                                                  unit_size=TREATMENT_SIZE,\n",
    "                                                  radius=RADIUS,\n",
    "                                                  pixel_value=pixel_value_iteration,\n",
    "                                                  constraint_layer=veg_constraint,\n",
    "                                                  # seed=seed # passing seed defined as diff number for every ee_treatment() run to randomize re-treatment groups\n",
    "                                                  ptOvershoot=PT_OVERSHOOT # we dial this in per scenario to get close to a 1:1 of Needed:Generated\n",
    "                                                  )\n",
    "    return treated_area_img,properties\n",
    "\n",
    "def remap_keys(d):\n",
    "  \"\"\"\n",
    "  solves a QA property problem for ee_treatments(): we have duplicate key names in the properties dict \n",
    "  when multiple ee_treatments() results are combined together.\n",
    "  So we create a unique identifier from 'yearInterval' and remove 'yearInterval' as a property itself.\n",
    "  example: 'PropertyName': PropertyValue becomes 'yearInterval1PropertyName': PropertyValue\n",
    "  \"\"\"\n",
    "  d = ee.Dictionary(d)\n",
    "  yearInterval = (ee.String('yearInterval')\n",
    "  .cat(ee.Number(ee.Dictionary(d).get('yearInterval')).format()))\n",
    "  keys = d.keys().slice(0,7) # depends on yearInterval property being in index position 7 of property list returned by ee_treatments()\n",
    "  values = d.values().slice(0,7)\n",
    "  remapped_keys = keys.map(lambda s: ee.String(yearInterval).cat(s))\n",
    "  return ee.Dictionary.fromLists(remapped_keys,values)\n",
    "\n",
    "def combine_dicts(d):\n",
    "   d = ee.Dictionary(d)\n",
    "   return d.keys().zip(d.values())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f0576d6f-47ed-445e-8b68-c94d8c382905",
   "metadata": {},
   "source": [
    "Create treated area ee.Images for every year interval on its assigned group of HUCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a2a9bc91",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "making new folder: projects/mas-gee/assets/treatment_scenarios/RunID18\n",
      "Export Started for projects/mas-gee/assets/treatment_scenarios/RunID18/NC_Fire_2m_20230707\n"
     ]
    }
   ],
   "source": [
    "# We use .map() to parallelize computation of treatment area images across the four 5-year intervals, \n",
    "# treating equal amount of area within each percentile ranking list, returning 4 treatment ee.Images, \n",
    "# then mosaicking them into one ee.Image for export\n",
    "\n",
    "# returns ee.List((ee.Image1,ee.Dictionary1),(ee.Image2,ee.Dictionary2),(etc))\n",
    "results = ee.List(zipped).map(treat_by_year_huc) \n",
    "\n",
    "# parse the images and properties out of the (image,properties) tuples by .map()ing over results\n",
    "images = results.map(lambda t: ee.List(t).get(0)) \n",
    "properties = results.map(lambda t: ee.List(t).get(1))\n",
    "\n",
    "# remap the keys of the properties so they are unique per ee_treatments() image run, then combine into one properties dict\n",
    "all_properties = properties.map(remap_keys).map(combine_dicts).flatten()\n",
    "# print(all_properties.getInfo())\n",
    "\n",
    "# mosaic all ee_treatment() images together and set the fixed properties \n",
    "output = ee.Image(ee.ImageCollection.fromImages(images).mosaic()).set(all_properties)\n",
    "\n",
    "# export\n",
    "desc = f'{RUNID}_{os.path.basename(output_image)}'\n",
    "\n",
    "folder = os.path.dirname(output_image)\n",
    "if check_exists(folder):\n",
    "    print(f'making new folder: {folder}')\n",
    "    os.popen(f'earthengine create folder {folder}').read()\n",
    "ee_treatments.export_img(output,desc,output_image,all_hucs.geometry(),30,'EPSG:5070')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.13 ('gee')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "5456adbb3bcda72c3f039b556ea69b5085f7c81f44c93a142b7acd907f0ad587"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}