example_forest_by_admin2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "59b88b8e-295e-4dee-bc19-cc7ae7d44740",
   "metadata": {},
   "source": [
    "# Example: Forested Area by Admin Level 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37c3de9a-1e5e-40d9-a016-784ee181e8b1",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "88a95b34-4d8e-4725-aaaf-227966e5cad1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "from earthengine_dask.core import ClusterGEE\n",
    "import google.auth\n",
    "import pandas as pd\n",
    "from pprint import pprint"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2019c415-4711-4299-a40c-4c5faafa6bf3",
   "metadata": {},
   "source": [
    "## Authenticate & Initialize Earth Engine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9026fd54-1285-4d65-a886-07115f452c8c",
   "metadata": {},
   "source": [
    "Get credentials and the GCP project ID, authenticating if necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7a15d6b4-2ff0-4a50-952d-4f2860de0624",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    credentials, project_id = google.auth.default()\n",
    "except google.auth.exceptions.DefaultCredentialsError:\n",
    "    !gcloud auth application-default login\n",
    "    credentials, project_id = google.auth.default()\n",
    "try:\n",
    "    ee.Initialize(credentials=credentials, project=project_id)\n",
    "except google.auth.exceptions.RefreshError:\n",
    "    !gcloud auth application-default login\n",
    "    credentials, project_id = google.auth.default()\n",
    "ee.Initialize(credentials=credentials, project=project_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84e59c6f-4d7f-4acf-ba20-6f15a43f98c4",
   "metadata": {},
   "source": [
    "# Input Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36265600-a446-4044-9efe-3422b5f0a21e",
   "metadata": {},
   "source": [
    "## Input: Forest Baseline\n",
    "\n",
    "This example will use the [European Commission Joint Research Centre's 2020 global map of forest cover](https://data.jrc.ec.europa.eu/dataset/10d1b337-b7d1-4938-a048-686c8185b290) for the forest baseline. The dataset is [available in Earth Engine](https://developers.google.com/earth-engine/datasets/catalog/JRC_GFC2020_V1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bc7df765-5c8b-45da-8bb6-705601e6743d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ic = ee.ImageCollection(\"JRC/GFC2020/V1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "439952f0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'bands': [],\n",
      " 'features': [{'bands': [{'crs': 'EPSG:4326',\n",
      "                          'crs_transform': [8.983152841195215e-05,\n",
      "                                            0,\n",
      "                                            -170.00005897568744,\n",
      "                                            0,\n",
      "                                            -8.983152841195215e-05,\n",
      "                                            80.03737653225383],\n",
      "                          'data_type': {'max': 255,\n",
      "                                        'min': 0,\n",
      "                                        'precision': 'int',\n",
      "                                        'type': 'PixelType'},\n",
      "                          'dimensions': [4007503, 1559941],\n",
      "                          'id': 'Map'}],\n",
      "               'id': 'JRC/GFC2020/V1/2020',\n",
      "               'properties': {'system:asset_size': 59186087403,\n",
      "                              'system:footprint': {'coordinates': [[-180, -90],\n",
      "                                                                   [180, -90],\n",
      "                                                                   [180, 90],\n",
      "                                                                   [-180, 90],\n",
      "                                                                   [-180, -90]],\n",
      "                                                   'type': 'LinearRing'},\n",
      "                              'system:index': '2020',\n",
      "                              'system:time_end': 1609372800000,\n",
      "                              'system:time_start': 1609372800000},\n",
      "               'type': 'Image',\n",
      "               'version': 1713492133769450.0}],\n",
      " 'id': 'JRC/GFC2020/V1',\n",
      " 'properties': {'max_mirrored_version': 1701018945816475},\n",
      " 'type': 'ImageCollection',\n",
      " 'version': 1713492133769450.0}\n"
     ]
    }
   ],
   "source": [
    "pprint(ic.getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f313f69c-04b5-46ca-b3b6-31458e3a4e74",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There is 1 image in the collection.\n"
     ]
    }
   ],
   "source": [
    "print(f'There is {ic.size().getInfo()} image in the collection.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2949b205-635e-4a48-9f33-cc8e1a8ec66f",
   "metadata": {},
   "source": [
    "... which we will use as the forest baseline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a26ba1fd-1116-49b3-a908-d9189696fd63",
   "metadata": {},
   "outputs": [],
   "source": [
    "jrc_forest_baseline = ic.first()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d85a36cf",
   "metadata": {},
   "source": [
    "Verify that the image is binary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "59ccc768",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Map_max': 1, 'Map_min': 0}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jrc_forest_baseline.unmask().reduceRegion(\n",
    "    reducer=ee.Reducer.minMax(),\n",
    "    geometry=ee.Geometry.BBox(10, 10, 11, 11),\n",
    "    maxPixels=1e10,\n",
    ").getInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d39a93e-db34-407d-80a2-134308b480f2",
   "metadata": {},
   "source": [
    "Looking at the projection information, the image is in decimal degrees of latitude and longitude (EPSG:4326)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6a41140f-bb9b-4c7e-909c-b8212ef50289",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'crs': 'EPSG:4326',\n",
      " 'transform': [8.983152841195215e-05,\n",
      "               0,\n",
      "               -170.00005897568744,\n",
      "               0,\n",
      "               -8.983152841195215e-05,\n",
      "               80.03737653225383],\n",
      " 'type': 'Projection'}\n"
     ]
    }
   ],
   "source": [
    "proj_info = jrc_forest_baseline.projection().getInfo()\n",
    "pprint(proj_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "97270963-97cb-4416-8fbd-e1ffb51e69af",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The nominal scale (at the equator) is 10 meters/pixel.\n"
     ]
    }
   ],
   "source": [
    "print(f'The nominal scale (at the equator) is '\n",
    "      f'{jrc_forest_baseline.projection().nominalScale().getInfo()} meters/pixel.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "930c34b3",
   "metadata": {},
   "source": [
    "## Input: Administrative Boundaries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c8fe28b-cd0b-4a33-aa8c-5df008e4a7ec",
   "metadata": {},
   "source": [
    "### GeoBoundaries ADM2\n",
    "\n",
    "We will use the municipal level (ADM2) boundaries provided by the [geoBoundaries](https://www.geoboundaries.org/) global database of political administrative boundaries v6.0, which is also [available in Earth Engine](https://developers.google.com/earth-engine/datasets/catalog/WM_geoLab_geoBoundaries_600_ADM2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c5073cf1-abdd-4413-a37e-99dbdd71a79f",
   "metadata": {},
   "outputs": [],
   "source": [
    "admin_fc = ee.FeatureCollection(\"WM/geoLab/geoBoundaries/600/ADM2\")\n",
    "# admin_fc = admin_fc.filter(ee.Filter.eq('shapeGroup', 'CAN'))\n",
    "# admin_fc = admin_fc.filter(ee.Filter.eq('shapeName', 'Colorado'))\n",
    "# admin_fc = admin_fc.filter(ee.Filter.eq('shapeName', 'Boulder'))\n",
    "\n",
    "# roi = ee.Geometry.Polygon(\n",
    "#         [[[-109.01952260759319, 40.971552045695994],\n",
    "#           [-109.01952260759319, 37.01127149086416],\n",
    "#           [-101.99925893571819, 37.01127149086416],\n",
    "#           [-101.99925893571819, 40.971552045695994]]], None, False)\n",
    "# admin_fc = admin_fc.filterBounds(roi)\n",
    "\n",
    "# admin_fc_props = None  # i.e. use all the properties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c34921bb-7257-4b49-b024-929e2f296427",
   "metadata": {},
   "source": [
    "There are quite a few features in the collection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "cdaa4d75-a6eb-4381-a916-9f4c88b02528",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 49617 features in the collection.\n"
     ]
    }
   ],
   "source": [
    "print(f'There are {admin_fc.size().getInfo()} features in the collection.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "521c96dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'shapeGroup': 'VAT',\n",
       " 'shapeID': '52476719B53484862947726',\n",
       " 'shapeName': 'Vatican City',\n",
       " 'shapeType': 'ADM0'}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "admin_fc.first().getInfo()['properties']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "53a0dda6-71fa-423d-9b21-f6ad1bb50017",
   "metadata": {},
   "outputs": [],
   "source": [
    "# admin.aggregate_histogram('shapeName').getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a9177880",
   "metadata": {},
   "outputs": [],
   "source": [
    "# geoboundaries_id_field = 'shapeID'\n",
    "# admin_fc_list = admin_fc.aggregate_histogram(geoboundaries_id_field).getInfo().keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cdf35fe-6240-46eb-bf00-e2950b997fd8",
   "metadata": {},
   "source": [
    "# Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bcd5ea2-76ca-424f-98eb-8bd46e88ade2",
   "metadata": {},
   "source": [
    "Define a function that calculates the forested area, and adds it back to the feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c520a58c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_area_base(\n",
    "        img,\n",
    "        feat,\n",
    "        properties=None,\n",
    "        max_pixels=1e10,\n",
    "    ):\n",
    "    stats_sum_temp = (\n",
    "        img.multiply(\n",
    "            ee.Image.pixelArea()\n",
    "        ).reduceRegion(\n",
    "            geometry=feat.geometry(),\n",
    "            reducer=ee.Reducer.sum(),\n",
    "            maxPixels=max_pixels,\n",
    "        )\n",
    "    )    \n",
    "    results = img.bandNames().map(\n",
    "        lambda band : ee.Number(stats_sum_temp.get(band)).round().multiply(1e-6)\n",
    "    ).getInfo()\n",
    "    # Return results for just the first band.\n",
    "    prop_dict = feat.toDictionary(properties).set('area_km2', results[0])\n",
    "    return ee.Dictionary(prop_dict).getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "83b0dbbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Aggreate a section at a time\n",
    "def get_area(\n",
    "        img,\n",
    "        feat,\n",
    "        properties=None,\n",
    "        max_pixels=1e10,\n",
    "        cover_grid_scale=16,\n",
    "    ):\n",
    "    feat = ee.Feature(feat)\n",
    "    \n",
    "    proj = img.projection().scale(2**cover_grid_scale, 2**cover_grid_scale)\n",
    "\n",
    "    covering_grid = feat.geometry().coveringGrid(proj)\n",
    "    # print('covering_grid.size() = ', covering_grid.size().getInfo())\n",
    "\n",
    "    grids = covering_grid.toList(1000)\n",
    "    geoms = grids.map(lambda f : ee.Feature(f).geometry())\n",
    "    area_total = 0\n",
    "    for grid_geom in geoms.getInfo():\n",
    "        f_grid = ee.Feature(grid_geom)\n",
    "        # TODO: test sensitivity to proj=img.projection(), maxError, reversing feat/f_grid\n",
    "        feature_piece = f_grid.intersection(feat, maxError=0.01, proj=img.projection())\n",
    "        piece_area = get_area_base(\n",
    "            img=img,\n",
    "            feat=feature_piece,\n",
    "            properties=properties,\n",
    "            max_pixels=max_pixels,\n",
    "        )\n",
    "        \n",
    "        # print('piece_area', piece_area.get('area_km2'))\n",
    "        area_total += piece_area.get('area_km2')\n",
    "    prop_dict = feat.toDictionary(properties).set('area_km2', area_total)\n",
    "    return ee.Dictionary(prop_dict).getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d17caf56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'area_km2': 12844.834378, 'shapeGroup': 'CAN', 'shapeID': '811477B35801915414510', 'shapeName': 'Muskoka--Kawarthas', 'shapeType': 'ADM2'}\n"
     ]
    }
   ],
   "source": [
    "def process_features(\n",
    "        img, \n",
    "        fc,\n",
    "        cover_grid_scale=16,\n",
    "    ):\n",
    "    id_list = fc.aggregate_array('system:index').getInfo()\n",
    "    for feature_id in id_list:\n",
    "        feat = fc.filter(ee.Filter.eq('system:index', feature_id)).first()\n",
    "        print(\n",
    "            get_area(\n",
    "                img,\n",
    "                feat,\n",
    "                properties=None,\n",
    "                cover_grid_scale=cover_grid_scale\n",
    "            )\n",
    "        )\n",
    "\n",
    "# process_features(img=jrc_forest_baseline, fc=admin_fc.filter(ee.Filter.eq('system:index', '00000000000000001c3a')))  # Icatu, Brazil\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(ee.Filter.eq('system:index', '00010000000000001fce'))\n",
    ")  # midsize Canada region\n",
    "# process_features(\n",
    "#     img=jrc_forest_baseline,\n",
    "#     fc=admin_fc.filter(\n",
    "#         ee.Filter.eq('system:index', '00010000000000001ff3')\n",
    "#     ),\n",
    "#     covering_grid_scale=16,\n",
    "# )  # Very large Northern Canada\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da604ce6",
   "metadata": {},
   "source": [
    "## Very large Northern Canada (00010000000000001ff3)\n",
    "\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00010000000000001ff3')\n",
    "    ),\n",
    "    covering_grid_scale=17,\n",
    ")  # Very large Northern Canada\n",
    "```\n",
    "Output (10 min)\n",
    "```\n",
    "covering_grid.size() =  8\n",
    "piece_area 1134.748005\n",
    "piece_area 7998.50299\n",
    "piece_area 2.614817\n",
    "piece_area 0\n",
    "piece_area 2723.368936\n",
    "piece_area 149812.669808\n",
    "piece_area 185034.036546\n",
    "piece_area 59548.806166999995\n",
    "{'area_km2': 406254.7472689999, 'shapeGroup': 'CAN', 'shapeID': '811477B98274245708554', 'shapeName': 'Northwest Territories / Terri*', 'shapeType': 'ADM2'}\n",
    "```\n",
    "\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00010000000000001ff3')\n",
    "    ),\n",
    "    covering_grid_scale=16,\n",
    ")  # Very large Northern Canada\n",
    "```\n",
    "Output (11 mins)\n",
    "```\n",
    "covering_grid.size() =  23\n",
    "piece_area 0\n",
    "piece_area 0.11370999999999999\n",
    "piece_area 0\n",
    "piece_area 0\n",
    "piece_area 1134.748005\n",
    "piece_area 4926.709285\n",
    "piece_area 3071.793705\n",
    "piece_area 2.459267\n",
    "piece_area 0.041839999999999995\n",
    "piece_area 0\n",
    "piece_area 2723.368936\n",
    "piece_area 37584.488913\n",
    "piece_area 80671.838106\n",
    "piece_area 57368.787766999994\n",
    "piece_area 32780.717534999996\n",
    "piece_area 4585.2317809999995\n",
    "piece_area 159.044887\n",
    "piece_area 225.228455\n",
    "piece_area 31331.114333999998\n",
    "piece_area 57098.651902\n",
    "piece_area 37785.879343\n",
    "piece_area 44006.422944\n",
    "piece_area 10798.106555999999\n",
    "{'area_km2': 406254.747271, 'shapeGroup': 'CAN', 'shapeID': '811477B98274245708554', 'shapeName': 'Northwest Territories / Terri*', 'shapeType': 'ADM2'}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3c8c1b7",
   "metadata": {},
   "source": [
    "## Hawaii 4 islands (00020000000000000e0b)\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00020000000000000e0b')\n",
    "    ),\n",
    "    covering_grid_scale=14,\n",
    ")  # Hawaii 4 islands\n",
    "```\n",
    "Output:\n",
    "```\n",
    "covering_grid.size() =  4\n",
    "piece_area 25.086164999999998\n",
    "piece_area 0.124369\n",
    "piece_area 246.526116\n",
    "piece_area 834.487941\n",
    "{'area_km2': 1106.224591, 'shapeGroup': 'USA', 'shapeID': '52423323B70973006487009', 'shapeName': 'Maui', 'shapeType': 'ADM2'}\n",
    "```\n",
    "\n",
    "\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00020000000000000e0b')\n",
    "    ),\n",
    "    covering_grid_scale=15,\n",
    ")  # Hawaii 4 islands\n",
    "```\n",
    "Output:\n",
    "```\n",
    "covering_grid.size() =  2\n",
    "piece_area 25.210534\n",
    "piece_area 1081.014056\n",
    "{'area_km2': 1106.22459, 'shapeGroup': 'USA', 'shapeID': '52423323B70973006487009', 'shapeName': 'Maui', 'shapeType': 'ADM2'}\n",
    "```\n",
    "\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00020000000000000e0b')\n",
    "    ),\n",
    "    covering_grid_scale=16,\n",
    ")  # Hawaii 4 islands\n",
    "```\n",
    "Output:\n",
    "```\n",
    "covering_grid.size() =  2\n",
    "piece_area 25.210534\n",
    "piece_area 1081.014056\n",
    "{'area_km2': 1106.22459, 'shapeGroup': 'USA', 'shapeID': '52423323B70973006487009', 'shapeName': 'Maui', 'shapeType': 'ADM2'}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49c49f55",
   "metadata": {},
   "source": [
    "## larger Canada region test\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00010000000000001fd2')\n",
    "    ),\n",
    "    covering_grid_scale=17,\n",
    ")  # larger Canada region\n",
    "```\n",
    "Output\n",
    "```\n",
    "covering_grid.size() =  1\n",
    "piece_area 222307.474921\n",
    "{'area_km2': 222307.474921, 'shapeGroup': 'CAN', 'shapeID': '811477B61209589162148', 'shapeName': 'Northeast / Nord-est', 'shapeType': 'ADM2'}\n",
    "```\n",
    "\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00010000000000001fd2')\n",
    "    ),\n",
    "    covering_grid_scale=16,\n",
    ")  # larger Canada region\n",
    "```\n",
    "Output\n",
    "```\n",
    "covering_grid.size() =  4\n",
    "piece_area 17097.961574\n",
    "piece_area 8281.468789\n",
    "piece_area 111576.56083799999\n",
    "piece_area 85351.483719\n",
    "{'area_km2': 222307.47492, 'shapeGroup': 'CAN', 'shapeID': '811477B61209589162148', 'shapeName': 'Northeast / Nord-est', 'shapeType': 'ADM2'}\n",
    "```\n",
    "\n",
    "```\n",
    "process_features(\n",
    "    img=jrc_forest_baseline,\n",
    "    fc=admin_fc.filter(\n",
    "        ee.Filter.eq('system:index', '00010000000000001fd2')\n",
    "    ),\n",
    "    covering_grid_scale=15,\n",
    ")  # larger Canada region\n",
    "```\n",
    "output\n",
    "```\n",
    "covering_grid.size() =  10\n",
    "piece_area 8160.096154999999\n",
    "piece_area 8937.865419\n",
    "piece_area 8281.468789\n",
    "piece_area 18003.483688\n",
    "piece_area 59734.079387\n",
    "piece_area 41074.933781\n",
    "piece_area 796.078213\n",
    "piece_area 33042.919551\n",
    "piece_area 38216.286307999995\n",
    "piece_area 6060.2636299999995\n",
    "{'area_km2': 222307.47492100002, 'shapeGroup': 'CAN', 'shapeID': '811477B61209589162148', 'shapeName': 'Northeast / Nord-est', 'shapeType': 'ADM2'}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c963133c",
   "metadata": {},
   "outputs": [],
   "source": [
    "## problematic polygon results\n",
    "\n",
    "## covering_grid_scale = 17 (? min)\n",
    "# covering_grid.size() =  17\n",
    "\n",
    "## covering_grid_scale = 16 (30 min)\n",
    "# covering_grid.size() =  44\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0.008088\n",
    "# piece_area 0.027437999999999997\n",
    "# piece_area 0.098663\n",
    "# piece_area 1.14146\n",
    "# piece_area 0.035799\n",
    "# piece_area 0.178333\n",
    "# piece_area 0.694444\n",
    "# piece_area 0.077383\n",
    "# piece_area 0.005212\n",
    "# piece_area 20.209909\n",
    "# piece_area 200.61459\n",
    "# piece_area 15.432127999999999\n",
    "# piece_area 46.374755\n",
    "# piece_area 20.339371\n",
    "# piece_area 0.6439969999999999\n",
    "# piece_area 0.9231809999999999\n",
    "# piece_area 1.8692369999999998\n",
    "# piece_area 4.3092\n",
    "# piece_area 0.568255\n",
    "# piece_area 0.015165999999999999\n",
    "# piece_area 3888.65222\n",
    "# piece_area 2585.778612\n",
    "# piece_area 0.123474\n",
    "# piece_area 0\n",
    "# piece_area 0.128498\n",
    "# piece_area 0.007783\n",
    "# piece_area 77.314599\n",
    "# piece_area 540.940532\n",
    "# {'area_km2': 7406.512327, 'shapeGroup': 'CAN', 'shapeID': '811477B26547766372806', 'shapeName': 'Nunavut', 'shapeType': 'ADM2'}\n",
    "\n",
    "## covering_grid_scale = 15 (60 min)\n",
    "# covering_grid.size() =  134\n",
    "# {'area_km2': 7406.512329000001, 'shapeGroup': 'CAN', 'shapeID': '811477B26547766372806', 'shapeName': 'Nunavut', 'shapeType': 'ADM2'}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee0134f8",
   "metadata": {},
   "source": [
    "\n",
    "## covering_grid_scale = 17 (? min)\n",
    "# covering_grid.size() =  17\n",
    "\n",
    "## covering_grid_scale = 16 (30 min)\n",
    "# covering_grid.size() =  44\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0\n",
    "# piece_area 0.008088\n",
    "# piece_area 0.027437999999999997\n",
    "# piece_area 0.098663\n",
    "# piece_area 1.14146\n",
    "# piece_area 0.035799\n",
    "# piece_area 0.178333\n",
    "# piece_area 0.694444\n",
    "# piece_area 0.077383\n",
    "# piece_area 0.005212\n",
    "# piece_area 20.209909\n",
    "# piece_area 200.61459\n",
    "# piece_area 15.432127999999999\n",
    "# piece_area 46.374755\n",
    "# piece_area 20.339371\n",
    "# piece_area 0.6439969999999999\n",
    "# piece_area 0.9231809999999999\n",
    "# piece_area 1.8692369999999998\n",
    "# piece_area 4.3092\n",
    "# piece_area 0.568255\n",
    "# piece_area 0.015165999999999999\n",
    "# piece_area 3888.65222\n",
    "# piece_area 2585.778612\n",
    "# piece_area 0.123474\n",
    "# piece_area 0\n",
    "# piece_area 0.128498\n",
    "# piece_area 0.007783\n",
    "# piece_area 77.314599\n",
    "# piece_area 540.940532\n",
    "# {'area_km2': 7406.512327, 'shapeGroup': 'CAN', 'shapeID': '811477B26547766372806', 'shapeName': 'Nunavut', 'shapeType': 'ADM2'}\n",
    "\n",
    "## covering_grid_scale = 15 (60 min)\n",
    "# covering_grid.size() =  134\n",
    "# {'area_km2': 7406.512329000001, 'shapeGroup': 'CAN', 'shapeID': '811477B26547766372806', 'shapeName': 'Nunavut', 'shapeType': 'ADM2'}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f35a3eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('med region in Canada')\n",
    "print(f'v2 {12844.829692}')\n",
    "print(f'v3 {9017.708524 + 3827.1211359999998}')\n",
    "print(f'v3 {9017.714211 + 3827.120167} proj=proj')\n",
    "print(f'v3 {9017.714211 + 3827.120167} proj=img.projection()')\n",
    "print('feature_piece = feat.intersection(f_grid, maxError=0.01, proj=img.projection())')\n",
    "print(f'v3 {9017.714071 + 3827.120203}')\n",
    "\n",
    "print('larger region in Canada')\n",
    "print(f'v2 {222307.47882699998}')\n",
    "print(f'v3 {17097.953794 + 8281.468525 + 111576.568655 + 85351.48794199999}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdc66b67-ad59-465a-b63e-ab3317b00f4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Use for debugging to count the total area, rather than forested area\n",
    "# forest_baseline = forest_baseline.unmask().multiply(0).add(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a377039f-4792-4d12-85a1-24c8b5d58b91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it out.\n",
    "# region_list = admin.aggregate_array('shapeID').distinct().getInfo()\n",
    "# region_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1ef0f78-9fc1-404c-9132-0203c690799f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # tileScale=1\n",
    "# get_area(img=forest_baseline, shape_id='42512837B26705409874577')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7690cfb9-2cac-447e-b235-45ad5782c413",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # tileScale=16\n",
    "# get_area(img=forest_baseline, shape_id='42512837B26705409874577')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50373e8c-cbe0-445d-873c-164edf14bc26",
   "metadata": {},
   "source": [
    "## Start Dask Cluster\n",
    "\n",
    "Start up a Earth Engine enabled cluster. This may take a few minutes to complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6693099f-82e7-4f0e-bd2a-a4da6dbcb91a",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = ClusterGEE(\n",
    "    name='test-cluster-forest-by-admin-temp4',\n",
    "    n_workers=5,\n",
    "    worker_cpu=8,\n",
    "    # spot_policy=\"spot_with_fallback\",\n",
    "    region='us-west1',\n",
    "    idle_timeout=\"4 hours\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a2c7acb-0564-4e3b-90ab-0e56134ab8cb",
   "metadata": {},
   "source": [
    "Retrieve a client for the cluster, and display it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02bb4842-09e3-4fb0-8899-97514ce283ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "client = cluster.get_client()\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "457b3b8c-cd92-4215-99e7-86fcdb1258a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problematic shape?\n",
    "#region_list = admin.aggregate_array('shapeID').distinct().getInfo()\n",
    "# region_list = ['42512837B26705409874577']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "248653a8-60fb-4b6c-87ec-054d1ec548c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create and submit jobs among the workers.\n",
    "# Allow for retries to handle \"Too many concurrent aggregations.\" errors\n",
    "submitted_jobs = [\n",
    "    client.submit(\n",
    "            get_area,  # function\n",
    "            binary_image, admin_fc, id_field, region_id, properties,  # function parameters\n",
    "            retries=1\n",
    "        )\n",
    "    for region_id in admin_regions\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72089446-90c5-4abd-8cc6-95cb9ad47dd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Debug issue with large regions by trying a problematic shape_id\n",
    "# submitted_jobs = [\n",
    "#     {\n",
    "#         'shape_id': '42512837B26705409874577',\n",
    "#         'tile_scale': tile_scale,\n",
    "#         'area':client.submit(\n",
    "#             get_area, forest_baseline, '42512837B26705409874577',\n",
    "#             retries=1\n",
    "#         )\n",
    "#     }\n",
    "# ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c68ab0e-b2cb-476e-a07f-865227f90e19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for job in submitted_jobs:\n",
    "#     if job.status in ['error']:\n",
    "#         future = job\n",
    "#         print(job)\n",
    "#         print(future.exception())\n",
    "\n",
    "jobs_with_errors = [\n",
    "    (job, job.exception(), )\n",
    "    for job in submitted_jobs\n",
    "    if job.status in ['error']\n",
    "]\n",
    "jobs_with_errors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dd1b104",
   "metadata": {},
   "source": [
    "## Get finished jobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc9c9611",
   "metadata": {},
   "outputs": [],
   "source": [
    "finished_jobs = [job for job in submitted_jobs if job.status=='finished']\n",
    "print(len(finished_jobs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b04c90ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "finished_results = client.gather(finished_jobs)\n",
    "finished_results[-5:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eac249a2-b763-414e-af35-063f916ca184",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(finished_results)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "396a1143",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show regions with the maximum alert area\n",
    "df.sort_values('area_km2').tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f62b01f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the results to a CSV file\n",
    "with open('output.csv', 'w') as out:\n",
    "    df.to_csv(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc5919d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#client.shutdown()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}