Unique_Zonal_Stats.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "9eb31b26-385e-498a-a19e-f44f150fcbb9",
      "metadata": {
        "id": "9eb31b26-385e-498a-a19e-f44f150fcbb9"
      },
      "source": [
        "# Purpose\n",
        "We are using GEE to calculate zonal statistics of EVI and LAI over the corn in California. We will masking them by a CDL corn layer but this could be a crop mask and exporting the data as a table.\n",
        "\n",
        "The unique zonal statistic for this task, using EVI as the example, is as follows:\n",
        "\n",
        "\\begin{align}\n",
        " \\frac{∑(EVI * Mask)}{∑(Mask)}\n",
        "    \\end{align}\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "365d8460-1946-40cf-9f58-c6bbcbea4f57",
      "metadata": {
        "id": "365d8460-1946-40cf-9f58-c6bbcbea4f57"
      },
      "source": [
        "Typical run me first."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3112c194-4d17-4f39-a1b6-224fea13bbb6",
      "metadata": {
        "tags": [],
        "id": "3112c194-4d17-4f39-a1b6-224fea13bbb6"
      },
      "outputs": [],
      "source": [
        "# Import ee\n",
        "import ee\n",
        "\n",
        "# Trigger the authentication flow.\n",
        "ee.Authenticate()\n",
        "\n",
        "# If tagging work\n",
        "# ee.data.setDefaultWorkloadTag('') # name the process in the cell here, for best results put in each cell"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5047e317-30c3-4fd5-9da3-98eb842ea265",
      "metadata": {
        "id": "5047e317-30c3-4fd5-9da3-98eb842ea265"
      },
      "source": [
        "# Imports\n",
        "Bring in some basic imports needed."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "179a02b6-bfc5-411e-990b-2c8eb2ff933b",
      "metadata": {
        "id": "179a02b6-bfc5-411e-990b-2c8eb2ff933b"
      },
      "source": [
        "## Geemap\n",
        "If you are wanting an interactive map, this needs to be done with every notebook and every time the notebook has sat offline for a while."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "711c01e6-47a4-40f1-b60e-09c1bc7155ae",
      "metadata": {
        "tags": [],
        "id": "711c01e6-47a4-40f1-b60e-09c1bc7155ae"
      },
      "outputs": [],
      "source": [
        "pip install geemap"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d3b3c75e-7d4a-45b9-b185-2671f5a3089d",
      "metadata": {
        "id": "d3b3c75e-7d4a-45b9-b185-2671f5a3089d"
      },
      "source": [
        "## Folium\n",
        "These maps are not truely interactive (i.e. can't inspect pixel values) but these are the most stable maps."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "266d370a-3cd4-4bfe-9778-7a091b46d629",
      "metadata": {
        "tags": [],
        "id": "266d370a-3cd4-4bfe-9778-7a091b46d629"
      },
      "outputs": [],
      "source": [
        "# Import folium for viewing maps\n",
        "import os\n",
        "import folium\n",
        "import geemap\n",
        "from datetime import datetime\n",
        "from IPython.display import Image\n",
        "#geemap.update_package()\n",
        "# Define the URL format used for Earth Engine generated map tiles.\n",
        "EE_TILES = 'https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}'\n",
        "\n",
        "# Initialize the library.\n",
        "ee.Initialize(project='') # put your gee project here!!!\n",
        "\n",
        "print('Folium version: ' + folium.__version__)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "448f81d6-232e-41c6-a799-771d7277b7dc",
      "metadata": {
        "id": "448f81d6-232e-41c6-a799-771d7277b7dc"
      },
      "source": [
        "## Displaying EE layers\n",
        "This function is similar to the one that Google presented in the tutorials but it adds the functionality of displaying image collections, geometries, and feature collections. This is if you are planning to use the folium map interface."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "64730e9c-1b4b-4c83-98d2-b5e6d7cbf803",
      "metadata": {
        "tags": [],
        "id": "64730e9c-1b4b-4c83-98d2-b5e6d7cbf803"
      },
      "outputs": [],
      "source": [
        "# Define a method for displaying Earth Engine image tiles on a folium map.\n",
        "def add_ee_layer(self, ee_object, vis_params, name):\n",
        "    try:\n",
        "        # Display ee.Image()\n",
        "        if isinstance(ee_object, ee.image.Image):\n",
        "            map_id_dict = ee.Image(ee_object).getMapId(vis_params)\n",
        "            folium.raster_layers.TileLayer(\n",
        "            tiles = map_id_dict['tile_fetcher'].url_format,\n",
        "            attr = 'Google Earth Engine',\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "            ).add_to(self)\n",
        "        # Display ee.ImageCollection()\n",
        "        elif isinstance(ee_object, ee.imagecollection.ImageCollection):\n",
        "            ee_object_new = ee_object.mosaic()\n",
        "            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n",
        "            folium.raster_layers.TileLayer(\n",
        "            tiles = map_id_dict['tile_fetcher'].url_format,\n",
        "            attr = 'Google Earth Engine',\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "            ).add_to(self)\n",
        "        # Display ee.Geometry()\n",
        "        elif isinstance(ee_object, ee.geometry.Geometry):\n",
        "            folium.GeoJson(\n",
        "            data = ee_object.getInfo(),\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "        ).add_to(self)\n",
        "        # Display ee.FeatureCollection()\n",
        "        elif isinstance(ee_object, ee.featurecollection.FeatureCollection):\n",
        "            ee_object_new = ee.Image().paint(ee_object, 0, 2)\n",
        "            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n",
        "            folium.raster_layers.TileLayer(\n",
        "            tiles = map_id_dict['tile_fetcher'].url_format,\n",
        "            attr = 'Google Earth Engine',\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "        ).add_to(self)\n",
        "\n",
        "    except:\n",
        "        print(\"Could not display {}\".format(name))\n",
        "\n",
        "# Add EE drawing method to folium.\n",
        "folium.Map.add_ee_layer = add_ee_layer\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5155e663-b336-4eda-a8e7-5d8f092a63c2",
      "metadata": {
        "id": "5155e663-b336-4eda-a8e7-5d8f092a63c2"
      },
      "source": [
        "# Date variables\n",
        "Here are the date variables that are going to be used in the experiment. If bringing in a new satellite source, enter the satellites start date here."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8f8475c6-5cc2-4ab9-829b-e64eca9cf731",
      "metadata": {
        "tags": [],
        "id": "8f8475c6-5cc2-4ab9-829b-e64eca9cf731"
      },
      "outputs": [],
      "source": [
        "# Satellite start dates\n",
        "mydSTART = '2002-07-04' # Aqua\n",
        "Today = datetime.today()\n",
        "\n",
        "yearStart = 2022 # 2003\n",
        "yearEnd   = 2022 # 2023\n",
        "\n",
        "dateStart = 121 # 121 - May 1\n",
        "dateEnd   = 160 # 288 - October 15\n",
        "\n",
        "format = 'DDD'\n",
        "\n",
        "dateStartD = ee.Date.parse(format, str(dateStart)).update(yearStart)\n",
        "dateEndD = ee.Date.parse(format, str(dateEnd)).update(yearEnd)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "73daa4f5-5a91-4f7e-9f8b-9337801d2edc",
      "metadata": {
        "id": "73daa4f5-5a91-4f7e-9f8b-9337801d2edc"
      },
      "source": [
        "# Admintrative layers\n",
        "Add in some proprietary or open source admin layers."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "5fde5320-e5bd-460f-96f0-51d18b77b516",
      "metadata": {
        "tags": [],
        "id": "5fde5320-e5bd-460f-96f0-51d18b77b516"
      },
      "outputs": [],
      "source": [
        "# Create a folium/gee map object.\n",
        "my_map = folium.Map(location=[34,265], zoom_start=5, height=1000) # Location is [lat,long]\n",
        "Mapa = geemap.Map(center=[34,265], zoom=5, height=1000)\n",
        "\n",
        "# Admin layer\n",
        "coun = ee.FeatureCollection(\"WM/geoLab/geoBoundaries/600/ADM0\").filter(ee.Filter.eq('shapeName','United States'))\n",
        "prov = ee.FeatureCollection(\"WM/geoLab/geoBoundaries/600/ADM1\").filter(ee.Filter.eq('shapeName','California'))\n",
        "dist = ee.FeatureCollection(\"WM/geoLab/geoBoundaries/600/ADM2\").filter(ee.Filter.eq('shapeName','Santa Barbara'))\n",
        "\n",
        "multiStateList=ee.List(['Iowa', 'Kansas', 'Texas'])\n",
        "def filterMulti(ids):\n",
        "    return ee.FeatureCollection(ids.map(lambda x: ee.FeatureCollection(\"WM/geoLab/geoBoundaries/600/ADM1\").filterMetadata('shapeName','equals', x))).flatten()\n",
        "multiReg = filterMulti(multiStateList)\n",
        "\n",
        "# Viewing the imagery\n",
        "#my_map.add_ee_layer(coun,{},'Country')\n",
        "#my_map.add_ee_layer(prov,{},'Province')\n",
        "#my_map.add_ee_layer(dist,{},'District')\n",
        "#my_map.add_ee_layer(multiReg,{},'Multiple Regions')\n",
        "#my_map.add_child(folium.LayerControl())\n",
        "#display(my_map)\n",
        "\n",
        "# Adding the imagery to the interactive panel\n",
        "Mapa.addLayer(coun,{},'Country')\n",
        "Mapa.addLayer(multiReg,{},'Multiple Regions')\n",
        "Mapa.addLayer(prov,{},'Province')\n",
        "Mapa.addLayer(dist,{},'District')\n",
        "Mapa.addLayerControl()\n",
        "Mapa"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# CDL Reference"
      ],
      "metadata": {
        "id": "5Cu_zdd4I8f3"
      },
      "id": "5Cu_zdd4I8f3"
    },
    {
      "cell_type": "code",
      "source": [
        "Mapb = geemap.Map(center=[34,265], zoom=5, height=1000)\n",
        "\n",
        "refYear = 2022\n",
        "\n",
        "cdlyear = str(refYear)\n",
        "\n",
        "# Select the CDL based on year\n",
        "cdl_20xx = ee.Image(ee.ImageCollection(\"USDA/NASS/CDL\").filter(ee.Filter.eq('system:index', cdlyear)).first()).select('cropland').clip(prov)\n",
        "#print(cdl_20xx)\n",
        "\n",
        "# Extract specific crops\n",
        "corn = cdl_20xx.eq(1).updateMask(cdl_20xx.eq(1)).toFloat().clip(prov)\n",
        "soy = cdl_20xx.eq(5).updateMask(cdl_20xx.eq(5)).toFloat()\n",
        "cotton = cdl_20xx.eq(2).updateMask(cdl_20xx.eq(2)).toFloat()\n",
        "rice = cdl_20xx.eq(3).updateMask(cdl_20xx.eq(3)).toFloat()\n",
        "sorghum = cdl_20xx.eq(4).updateMask(cdl_20xx.eq(4)).toFloat()\n",
        "dwheat = cdl_20xx.eq(22).updateMask(cdl_20xx.eq(22)).toFloat()\n",
        "swheat = cdl_20xx.eq(23).updateMask(cdl_20xx.eq(23)).toFloat()\n",
        "wwheat = cdl_20xx.eq(24).updateMask(cdl_20xx.eq(24)).toFloat()\n",
        "sugarc = cdl_20xx.eq(25).updateMask(cdl_20xx.eq(25)).toFloat()\n",
        "\n",
        "# Add to map\n",
        "#Mapb.addLayer(cdl_20xx,{},'CDL')\n",
        "Mapb.addLayer(corn,{'palette':'ffd400'},'Corn')\n",
        "Mapb.addLayer(soy,{'palette':'267000'},'Soy')\n",
        "Mapb.addLayerControl()\n",
        "Mapb"
      ],
      "metadata": {
        "id": "cfPjGLjRI6kz"
      },
      "id": "cfPjGLjRI6kz",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "id": "abd9473a-0ce2-46ab-8a62-68db688c41af",
      "metadata": {
        "id": "abd9473a-0ce2-46ab-8a62-68db688c41af"
      },
      "source": [
        "# Load satellite imagery\n",
        "Here we load AQUA EVI and LAI but any environmental source can be loaded."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3563d1b4-4125-4c05-b992-bc95478d1315",
      "metadata": {
        "tags": [],
        "id": "3563d1b4-4125-4c05-b992-bc95478d1315"
      },
      "outputs": [],
      "source": [
        "# Load MYD13A1.061: Aqua Vegetation Indices 16-Day Global 500m\n",
        "myd13 = ee.ImageCollection(\"MODIS/061/MYD13A1\").filterDate(mydSTART, Today)\n",
        "\n",
        "# Load MYD15A2H.061: Aqua Leaf Area Index/FPAR 8-Day Global 500m\n",
        "myd15 = ee.ImageCollection(\"MODIS/061/MYD15A2H\").filterDate(mydSTART, Today)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a75b9ee3-7226-4e58-83d2-76bb3ac6e7c7",
      "metadata": {
        "id": "a75b9ee3-7226-4e58-83d2-76bb3ac6e7c7"
      },
      "source": [
        "# View/process satellite imagery\n",
        "This section is dedicated to visualizing the imagery as well as creating the variables needed for the zonal statistics section. If you are going to bring in a new source just copy and paste one of the code blocks and change the variables as appropriate."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8b871e9d-99f2-4c95-955c-e4afb618466c",
      "metadata": {
        "id": "8b871e9d-99f2-4c95-955c-e4afb618466c"
      },
      "source": [
        "## View Aqua MODIS EVI"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "5c5875f3-829a-4a44-981d-233fb5aae593",
      "metadata": {
        "tags": [],
        "id": "5c5875f3-829a-4a44-981d-233fb5aae593"
      },
      "outputs": [],
      "source": [
        "# Create a folium/gee map object.\n",
        "my_map = folium.Map(location=[34,265], zoom_start=5, height=1000) # Location is [lat,long]\n",
        "Map8a = geemap.Map(center=[34,265], zoom=5, height=1000)\n",
        "\n",
        "############Aqua MODIS MYD13 EVI############\n",
        "def findEvi(i):\n",
        "  iveg = i.select('EVI').multiply(0.0001)\n",
        "  return iveg\n",
        "\n",
        "eviFun = myd13.filterDate(dateStartD, dateEndD).map(findEvi)\n",
        "eviFunMed = ee.ImageCollection(eviFun).median()\n",
        "\n",
        "# Viewing the imagery\n",
        "#my_map.add_ee_layer(eviFunMed.clip(multiReg),{'palette':['ff0000','ffff00','00ff00']},'EVI median')\n",
        "#my_map.add_ee_layer(eviFun,{'min':0.19,'max':0.32},'EVI raw')\n",
        "#my_map.add_child(folium.LayerControl())\n",
        "#display(my_map)\n",
        "\n",
        "# Adding the imagery to the interactive panel\n",
        "Map8a.addLayer(eviFunMed.clip(prov),{'min':0,'max':0.9,'palette':['ff0000','ffff00','00ff00']},'EVI median')\n",
        "#Map8a.addLayer(eviFun,{'min':0.19,'max':0.32},'EVI raw')\n",
        "Map8a.addLayerControl()\n",
        "Map8a"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d86c7243-79b7-425d-99ce-2ad471a2edaf",
      "metadata": {
        "id": "d86c7243-79b7-425d-99ce-2ad471a2edaf"
      },
      "source": [
        "### Calculate EVI on Crop Mask"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f82fa47b-3a85-4d3c-9a73-8bf4380981e1",
      "metadata": {
        "id": "f82fa47b-3a85-4d3c-9a73-8bf4380981e1"
      },
      "outputs": [],
      "source": [
        "# Create a folium/gee map object.\n",
        "my_map = folium.Map(location=[34,265], zoom_start=5, height=1000) # location is [lat,long]\n",
        "Map8aa = geemap.Map(center=[34,265], zoom=5, height=1000)\n",
        "\n",
        "############Crop mask calculations############\n",
        "mask = ee.Image(corn) # Crop mask here or CDL crop layer\n",
        "maskFloat = mask.divide(100) # convert crop mask to float\n",
        "\n",
        "def findEviOnMask(i):\n",
        "  evi = i.select('EVI').multiply(0.0001) # extract the correct layer and apply the scaling factor\n",
        "  eviOnMa = evi.multiply(maskFloat)\n",
        "\n",
        "  joinFilter = ee.Filter.equals(\n",
        "  leftField='shapeID',\n",
        "  rightField='shapeID'\n",
        "  )\n",
        "\n",
        "  # Get summed EVI * floating mask per region per date\n",
        "  summedEOnM = eviOnMa.reduceRegions(\n",
        "      collection=prov, # add the shapefile here\n",
        "      reducer=ee.Reducer.sum(),\n",
        "      scale=100\n",
        "  ).select(**{\n",
        "    'propertySelectors': ['shapeGroup','shapeName','shapeID','sum'],\n",
        "    'newProperties': ['shapeGroup','shapeName','shapeID','sumEOnM']\n",
        "    })\n",
        "\n",
        "  # Get summed floating mask per region per date\n",
        "  summedMask = maskFloat.reduceRegions(\n",
        "      collection=prov, # add the shapefile here\n",
        "      reducer=ee.Reducer.sum(),\n",
        "      scale=100\n",
        "  ).select(**{\n",
        "    'propertySelectors': ['shapeID','sum'],\n",
        "    'newProperties': ['shapeID','sumMask']\n",
        "    })\n",
        "\n",
        "  joinHere = ee.Join.inner()\n",
        "  joinedJoin = joinHere.apply(summedEOnM,summedMask,joinFilter)\n",
        "  def calcDone(i):\n",
        "      getProp = ee.Feature(i.get('primary')).copyProperties(i.get('secondary'))\n",
        "      return getProp.set({'calcedVal': ee.Number(getProp.get('sumEOnM')).divide(ee.Number(getProp.get('sumMask')))}) # divide the sums\n",
        "  calculatedDistricts = joinedJoin.map(calcDone)\n",
        "\n",
        "  return calculatedDistricts\n",
        "\n",
        "eviOnMaFun = myd13.filterDate(dateStartD, dateEndD).map(findEviOnMask).flatten()\n",
        "eviOnMaFunMed = ee.ImageCollection(eviOnMaFun).median()\n",
        "\n",
        "# Viewing the imagery\n",
        "#my_map.add_ee_layer(eviFunMed.clip(multiReg),{'palette':['ff0000','ffff00','00ff00']},'EVI median')\n",
        "#my_map.add_ee_layer(eviFun,{'min':0.19,'max':0.32},'EVI raw')\n",
        "#my_map.add_child(folium.LayerControl())\n",
        "#display(my_map)\n",
        "\n",
        "# Adding the imagery to the interactive panel\n",
        "#Map8aa.addLayer(mask.clip(multiReg),{'min':0,'max':100,},'Crop mask')\n",
        "#Map8aa.addLayer(maskFloat.clip(multiReg),{'min':0,'max':1,},'Crop mask Float')\n",
        "#Map8aa.addLayer(eviFunMed.clip(multiReg),{'min':0,'max':1},'EVI Median')\n",
        "Map8aa.addLayer(eviOnMaFun,{},'EVI on Mask')\n",
        "Map8aa.addLayerControl()\n",
        "Map8aa\n",
        "\n",
        "#print(eviOnMaFun.getInfo())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3926ceec-bb61-4181-9080-d7af2aaaf3ca",
      "metadata": {
        "id": "3926ceec-bb61-4181-9080-d7af2aaaf3ca"
      },
      "source": [
        "#### Export EVI on Crop Mask"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d0dba5e9-2e82-4cfe-8f5e-1222a18da02f",
      "metadata": {
        "id": "d0dba5e9-2e82-4cfe-8f5e-1222a18da02f"
      },
      "outputs": [],
      "source": [
        "# Export the statistics\n",
        "taskExE = ee.batch.Export.table.toDrive(**{\n",
        "    'collection':eviOnMaFun,\n",
        "    'folder':'Zonal_Stats', # you need to put the name of a folder from your Google Drive here\n",
        "    'description':'EVIonCorn_season_2022', # name of the exported file\n",
        "    'fileFormat':'CSV',\n",
        "    'selectors':['system:index','shapeGroup','shapeName','shapeID','sumEOnM','sumMask','calcedVal'] # columns\n",
        "})\n",
        "\n",
        "taskExE.start()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "294f3e28-0a58-4e53-8642-90b7ee60e4c7",
      "metadata": {
        "id": "294f3e28-0a58-4e53-8642-90b7ee60e4c7"
      },
      "outputs": [],
      "source": [
        "# Check the task status\n",
        "task_status = taskExE.status()\n",
        "\n",
        "# If the task is running\n",
        "if task_status[\"state\"] == \"RUNNING\":\n",
        "  start_time = task_status[\"start_timestamp_ms\"]/(60000)\n",
        "  update_time = task_status['update_timestamp_ms']/(60000)\n",
        "  elapsed_time = update_time-(start_time)\n",
        "  hours, minutes = divmod(elapsed_time, 60)\n",
        "  # Print a message saying that the task is running and how long it has been running\n",
        "  print(\"The task is running. It has been running for {} hours and {} minutes.\".format(round(hours,1), round(minutes,1)))\n",
        "\n",
        "# If the task has finished\n",
        "elif task_status[\"state\"] == \"COMPLETED\":\n",
        "  # Print a message saying that the task has finished\n",
        "  print(\"The task has finished.\")\n",
        "\n",
        "# If the task is in an unknown state\n",
        "else:\n",
        "  # Print a message saying that the task is in an unknown state\n",
        "  print(\"The task is in an unknown state.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "00909f04-7e2e-47ca-b336-c80f3d27b1e9",
      "metadata": {
        "id": "00909f04-7e2e-47ca-b336-c80f3d27b1e9"
      },
      "source": [
        "## View Aqua MODIS LAI"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "37df022c-1424-46e2-baee-03504d4331ab",
      "metadata": {
        "id": "37df022c-1424-46e2-baee-03504d4331ab"
      },
      "outputs": [],
      "source": [
        "# Create a folium/gee map object.\n",
        "my_map = folium.Map(location=[34,265], zoom_start=5, height=1000) # Location is [lat,long]\n",
        "Map8b = geemap.Map(center=[34,265], zoom=5, height=1000)\n",
        "\n",
        "############Aqua MODIS MYD15A2H LAI############\n",
        "def findLai(i):\n",
        "  iveg = i.select('Lai_500m').multiply(0.1)\n",
        "  return iveg\n",
        "\n",
        "laiFun = myd15.filterDate(dateStartD, dateEndD).map(findLai)\n",
        "laiFunMed = ee.ImageCollection(laiFun).median()\n",
        "\n",
        "# Viewing the imagery\n",
        "#my_map.add_ee_layer(laiFunMed.clip(multiReg),{'palette':['ff0000','ffff00','00ff00']},'LAI median')\n",
        "#my_map.add_ee_layer(laiFun,{'min':0.19,'max':0.32},'LAI raw')\n",
        "#my_map.add_child(folium.LayerControl())\n",
        "#display(my_map)\n",
        "\n",
        "# Adding the imagery to the interactive panel\n",
        "Map8b.addLayer(laiFunMed.clip(prov),{'min':0,'max':3,'palette':['ff0000','ffff00','00ff00']},'LAI median')\n",
        "#Map8b.addLayer(laiFun,{'min':0.19,'max':0.32},'LAI raw')\n",
        "Map8b.addLayerControl()\n",
        "Map8b"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0f5ba2ac-61c0-4a47-b131-c498e56fac11",
      "metadata": {
        "id": "0f5ba2ac-61c0-4a47-b131-c498e56fac11"
      },
      "source": [
        "### Calculate LAI on Crop Mask"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "03837623-7569-4d48-babb-24f690903759",
      "metadata": {
        "id": "03837623-7569-4d48-babb-24f690903759"
      },
      "outputs": [],
      "source": [
        "# Create a folium/gee map object.\n",
        "my_map = folium.Map(location=[34,265], zoom_start=5, height=1000) # location is [lat,long]\n",
        "Map8aa = geemap.Map(center=[34,265], zoom=5, height=1000)\n",
        "\n",
        "############Crop mask calculations############\n",
        "mask = ee.Image(corn) # Crop mask here or CDL crop layer\n",
        "maskFloat = mask.divide(100) # convert crop mask to float\n",
        "\n",
        "def findLaiOnMask(i):\n",
        "  lai = i.select('Lai_500m').multiply(0.1) # extract the correct layer and apply the scaling factor\n",
        "  laiOnMa = lai.multiply(maskFloat)\n",
        "\n",
        "  joinFilter = ee.Filter.equals(\n",
        "  leftField='shapeID',\n",
        "  rightField='shapeID'\n",
        "  )\n",
        "\n",
        "  # Get summed LAI * floating mask per region per date\n",
        "  summedLOnM = laiOnMa.reduceRegions(\n",
        "      collection=prov, # add the shapefile here\n",
        "      reducer=ee.Reducer.sum(),\n",
        "      scale=100\n",
        "  ).select(**{\n",
        "    'propertySelectors': ['shapeGroup','shapeName','shapeID','sum'],\n",
        "    'newProperties': ['shapeGroup','shapeName','shapeID','sumLOnM']\n",
        "    })\n",
        "\n",
        "  # Get summed floating mask per region per date\n",
        "  summedMask = maskFloat.reduceRegions(\n",
        "      collection=prov, # add the shapefile here\n",
        "      reducer=ee.Reducer.sum(),\n",
        "      scale=100\n",
        "  ).select(**{\n",
        "    'propertySelectors': ['shapeID','sum'],\n",
        "    'newProperties': ['shapeID','sumMask']\n",
        "    })\n",
        "\n",
        "  # Applying the join and perform the calculations over each region\n",
        "  joinHere = ee.Join.inner()\n",
        "  joinedJoin = joinHere.apply(summedLOnM,summedMask,joinFilter)\n",
        "  def calcDone(i):\n",
        "      getProp = ee.Feature(i.get('primary')).copyProperties(i.get('secondary'))\n",
        "      return getProp.set({'calcedVal': ee.Number(getProp.get('sumLOnM')).divide(ee.Number(getProp.get('sumMask')))}) # divide the sums\n",
        "  calculatedDistricts = joinedJoin.map(calcDone)\n",
        "\n",
        "  return calculatedDistricts\n",
        "\n",
        "laiOnMaFun = myd15.filterDate(dateStartD, dateEndD).map(findLaiOnMask).flatten()\n",
        "laiOnMaFunMed = ee.ImageCollection(laiOnMaFun).median()\n",
        "\n",
        "# Viewing the imagery\n",
        "#my_map.add_ee_layer(laiFunMed.clip(multiReg),{'palette':['ff0000','ffff00','00ff00']},'LAI median')\n",
        "#my_map.add_ee_layer(laiFun,{'min':0.19,'max':0.32},'LAI raw')\n",
        "#my_map.add_child(folium.LayerControl())\n",
        "#display(my_map)\n",
        "\n",
        "# Adding the imagery to the interactive panel\n",
        "#Map8aa.addLayer(mask.clip(multiReg),{'min':0,'max':100,},'Crop mask')\n",
        "#Map8aa.addLayer(maskFloat.clip(multiReg),{'min':0,'max':1,},'Crop mask Float')\n",
        "#Map8aa.addLayer(laiFunMed.clip(multiReg),{'min':0,'max':1},'LAI Median')\n",
        "Map8aa.addLayer(laiOnMaFun,{},'LAI on Mask')\n",
        "Map8aa.addLayerControl()\n",
        "Map8aa\n",
        "\n",
        "#print(eviOnMaFun.getInfo())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a5e4e0aa-f325-4b28-887f-3f3d013d665b",
      "metadata": {
        "id": "a5e4e0aa-f325-4b28-887f-3f3d013d665b"
      },
      "source": [
        "#### Export LAI on Crop Mask"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "26528c49-08a4-457f-a6dd-6e5a314e1f2b",
      "metadata": {
        "id": "26528c49-08a4-457f-a6dd-6e5a314e1f2b"
      },
      "outputs": [],
      "source": [
        "# Export the statistics\n",
        "taskExL = ee.batch.Export.table.toDrive(**{\n",
        "    'collection':laiOnMaFun,\n",
        "    'folder':'Zonal_Stats', # you need to put the name of a folder from your Google Drive here\n",
        "    'description':'LAIonCorn_season_202', # name of the exported file\n",
        "    'fileFormat':'CSV',\n",
        "    'selectors':['system:index','shapeGroup','shapeName','shapeID','sumLOnM','sumMask','calcedVal'] # columns\n",
        "})\n",
        "\n",
        "taskExL.start()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "19dac59e-e874-4c03-bee9-f4274e4b8bdb",
      "metadata": {
        "id": "19dac59e-e874-4c03-bee9-f4274e4b8bdb"
      },
      "outputs": [],
      "source": [
        "# Check the task status\n",
        "task_status = taskExL.status()\n",
        "\n",
        "# If the task is running\n",
        "if task_status[\"state\"] == \"RUNNING\":\n",
        "  start_time = task_status[\"start_timestamp_ms\"]/(60000)\n",
        "  update_time = task_status['update_timestamp_ms']/(60000)\n",
        "  elapsed_time = update_time-(start_time)\n",
        "  hours, minutes = divmod(elapsed_time, 60)\n",
        "  # Print a message saying that the task is running and how long it has been running\n",
        "  print(\"The task is running. It has been running for {} hours and {} minutes.\".format(round(hours,1), round(minutes,1)))\n",
        "\n",
        "# If the task has finished\n",
        "elif task_status[\"state\"] == \"COMPLETED\":\n",
        "  # Print a message saying that the task has finished\n",
        "  print(\"The task has finished.\")\n",
        "\n",
        "# If the task is in an unknown state\n",
        "else:\n",
        "  # Print a message saying that the task is in an unknown state\n",
        "  print(\"The task is in an unknown state.\")"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "gee",
      "language": "python",
      "name": "gee"
    },
    "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.12.0"
    },
    "colab": {
      "provenance": [],
      "toc_visible": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}