1. Deforestation attribution (GEE)-Spatial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "05ceeb3c-2ef4-4cdb-b443-b8f90e9582e2",
   "metadata": {},
   "source": [
    "# Authentiacate GEE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "95e5b55d-3a4e-4ef4-bf53-cfd75fe6be7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This will track the datasets version if any change is made to the code\n",
    "GEE_version = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b46bf3c1-556c-47a8-8f58-9ba4251fd172",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<p>To authorize access needed by Earth Engine, open the following\n",
       "        URL in a web browser and follow the instructions:</p>\n",
       "        <p><a href=https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=64frUDoZxdZYcEEndhUl0DYYsbXas85LVktx7nTRlyU&tc=pL5wQ80Idu5U1FHAqMKyXHRpJHNSrc_SJzZAu7n50VA&cc=mYDrD-YDGEy_W-HUHnnJ3q0-YyLQ0ZEh02SCDl-Q70k>https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=64frUDoZxdZYcEEndhUl0DYYsbXas85LVktx7nTRlyU&tc=pL5wQ80Idu5U1FHAqMKyXHRpJHNSrc_SJzZAu7n50VA&cc=mYDrD-YDGEy_W-HUHnnJ3q0-YyLQ0ZEh02SCDl-Q70k</a></p>\n",
       "        <p>The authorization workflow will generate a code, which you should paste in the box below.</p>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Enter verification code:  4/1AQlEd8ymL_JdzJqlRXhDJNDKM_HUhlPZtFniEYmpw-V8mi23Tqm4yMMtScM\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Successfully saved authorization token.\n"
     ]
    }
   ],
   "source": [
    "import ee\n",
    "# Trigger the authentication flow for Google Earth Engine (GEE).\n",
    "ee.Authenticate()\n",
    "# Initialize the library.\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "229b7218-5ff1-4fd5-8077-3edf5b67725d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading all the nessesary packages\n",
    "import datetime\n",
    "import geetools\n",
    "import geemap\n",
    "import concurrent.futures\n",
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5163de8f-eb98-4d2b-8022-cf7bbd11bbb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(253,)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the feature collection (i.e., boundary shapefile in csv format) needed to simulate aggregation\n",
    "# Note: that this csv file will just contain countries names, but the actual boundary (be it first or second level admin boundary) that is being simulated should be uploaded in GEE\n",
    "GADM_countries = pd.read_csv(\"/home/chandra/backup/Chalmers/GADM_Boundaries.csv\")\n",
    "\n",
    "# Groupby and count the attributes/shapes in the dataframe, just to get a sense of number of distinct sub-boundaries within a country\n",
    "counts = GADM_countries.groupby('COUNTRY').size().reset_index(name='count')\n",
    "\n",
    "# Resulting dataframe with distinct Countries and unique sub-boundaries within them\n",
    "Countries = pd.DataFrame(counts[['COUNTRY', 'count']]).sort_values(by='count')\n",
    "\n",
    "# Check total countries that will be simulated\n",
    "GADM_countries['COUNTRY'].unique().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "39a9ea0d-0129-415e-9a07-7a3a335b643f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def GEE_loop(Admin_boundary, simulation = 'All'):\n",
    "    global geometry, uniqueValues_ADM2CODE\n",
    "    \n",
    "    ####################\n",
    "    # 1. Intial run values \n",
    "    ####################\n",
    "    \n",
    "    # Add end year , i.e., the year utpto which the attribution simulation runs\n",
    "    end_year = 2022\n",
    "    # Assuming tree cover greater than and equal to (≥) 25% is considered forest in Hansen dataset\n",
    "    Forest_threshold = 25\n",
    "    # Forestloss attribution to plantation, by assuming that pixels with startyear >2000 are due to deforestation after year 2000, the rest is considered under deforestation pre-2000's\n",
    "    Plantation_threshold_year = 2000\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    ####################\n",
    "    # 2. Call/Check latest datasets \n",
    "    ####################\n",
    "    \n",
    "    # Load the feature collection (i.e., shapefile) needed to simulate aggregation \n",
    "    GADM_admin_boundary_version = \"projects/lu-chandrakant/assets/GADM_dissolved\"   # This file is downloaded from GADM (version 4.1; https://gadm.org/data.html) and uploaded to google earth engine\n",
    "    \n",
    "    # Other raster datasets\n",
    "    Hansen_tree_cover_dataset_version = \"UMD/hansen/global_forest_change_2022_v1_10\" # Check the latest version of Hansen et al. 2013 (https://www.science.org/doi/10.1126/science.1244693p) at https://glad.earthengine.app/view/global-forest-change\n",
    "    Plantation_dataset_version = \"projects/lu-chandrakant/assets/Plantation_new/\"    # Check latest plantation dataset from Du et al. 2022 (https://www.nature.com/articles/s41597-022-01260-2) and upload to GEE\n",
    "    \n",
    "    # Update/Check MapBiomas dataset colllection (Integrated maps)\n",
    "    MapBiomas_dataset_version = {'Amazon': 'projects/mapbiomas-raisg/public/collection5/mapbiomas_raisg_panamazonia_collection5_integration_v1',                  #MapBiomas-Amazonia (1985-2022) \n",
    "                                  'Argentina': 'projects/mapbiomas-public/assets/argentina/collection1/mapbiomas_argentina_collection1_integration_v1',           #MapBiomas-Argentina (1998-2022)\n",
    "                                  'Atlantic forest': 'projects/mapbiomas_af_trinacional/public/collection3/mapbiomas_atlantic_forest_collection30_integration_v1',#MapBiomas-Atlantic forest (1985-2022)\n",
    "                                  'Brazil': 'projects/mapbiomas-workspace/public/collection8/mapbiomas_collection80_integration_v1',                              #MapBiomas-Brazil (1985-2022)\n",
    "                                  'Chaco': 'projects/mapbiomas-chaco/public/collection4/mapbiomas_chaco_collection4_integration_v1',                              #MapBiomas-Chaco (1985-2022)\n",
    "                                  'Chile': 'projects/mapbiomas-public/assets/chile/collection1/mapbiomas_chile_collection1_integration_v1',                       #MapBiomas-Chile (2000-2022)\n",
    "                                  'Colombia': 'projects/mapbiomas-public/assets/colombia/collection1/mapbiomas_colombia_collection1_integration_v1',              #MapBiomas-Colombia (1985-2022)\n",
    "                                  'Ecaudor': 'projects/mapbiomas-public/assets/ecuador/collection1/mapbiomas_ecuador_collection1_integration_v1',                 #MapBiomas-Ecaudor (1985-2022)\n",
    "                                  'Indonesia': 'projects/mapbiomas-indonesia/public/collection2/mapbiomas_indonesia_collection2_integration_v1',                  #MapBiomas-Indonesia (2000-2022)\n",
    "                                  'Pampa': 'projects/MapBiomas_Pampa/public/collection3/mapbiomas_pampa_collection3_integration_v1',                              #MapBiomas-Pampa (1985-2022)\n",
    "                                  'Paraguay': 'projects/mapbiomas-public/assets/paraguay/collection1/mapbiomas_paraguay_collection1_integration_v1',              #MapBiomas-Paraguay (1985-2022)\n",
    "                                  'Peru': 'projects/mapbiomas-public/assets/peru/collection2/mapbiomas_peru_collection2_integration_v1',                          #MapBiomas-Peru (1985-2022)\n",
    "                                  'Uruguay': 'projects/MapBiomas_Pampa/public/collection3/mapbiomas_uruguay_collection1_integration_v1',                          #MapBiomas-Uruguay (1985-2022)\n",
    "                                  'Venezuela': 'projects/mapbiomas-public/assets/venezuela/collection1/mapbiomas_venezuela_collection1_integration_v1',           #MapBiomas-Venezuela (1985-2022)\n",
    "                                  'Bolivia': 'projects/mapbiomas-public/assets/bolivia/collection1/mapbiomas_bolivia_collection1_integration_v1',                 #MapBiomas-Bolivia (1985-2021)\n",
    "                            }                                                                                                                  \n",
    "                                 \n",
    "    Soy_dataset_version = 'projects/glad/soy_annual_SA' # Update with latest version Soybean dataset from Song et al. 2020 (https://www.nature.com/articles/s41893-021-00729-z). \n",
    "                                                        # Presently contains Soy data till 2022. Data for 2021 is uploaded manually. Sometime the latest year data may have to be inputted manually (Code's inactive part in \"\"\"   \"\"\")\n",
    "\n",
    "    Cocoa_dataset_version = 'projects/ee-nk-cocoa/assets/cocoa_map_threshold_065'  # Latest dataset from Kalischek et al. 2023 (https://doi.org/10.1038/s43016-023-00751-8)\n",
    "    Oilpalm_Indonesia_version = 'projects/lu-chandrakant/assets/OilPalm_Indonesia' # Latest dataset from Gaveau et al. 2022 (https://doi.org/10.1371/journal.pone.0266178)\n",
    "    Oilpalm_Malaysia_version = 'projects/lu-chandrakant/assets/OilPalm_Indonesia_Malaysia'    # Latest dataset from Xu et al. 2020 (https://doi.org/10.5194/essd-12-847-2020)\n",
    "    Oilpalm_Global_version = \"BIOPAMA/GlobalOilPalm/v1\"                            # Latest dataset from Descals et al. 2021 (https://doi.org/10.5194/essd-13-1211-2021)\n",
    "    Coconut_dataset_version = \"projects/lu-chandrakant/assets/Coconut\"             # Latest dataset from Descals et al. 2023 (https://doi.org/10.5194/essd-15-3991-2023)\n",
    "    Sugarcane_dataset_version = \"projects/lu-chandrakant/assets/Sugarcane\"         # Latest dataset from Zheng et al. 2022 (https://doi.org/10.5194/essd-14-2065-2022)\n",
    "    Rapeseed_dataset_version = \"projects/lu-chandrakant/assets/Rapeseed\"           # Latest dataset from Han et al. 2021 (https://doi.org/10.5194/essd-13-2857-2021)\n",
    "    Rice_dataset_version = \"projects/lu-chandrakant/assets/Paddy_rice\"             # Latest dataset from Han et al. 2021 (https://doi.org/10.5194/essd-13-5969-2021)\n",
    "    Maize_China_dataset_version = \"projects/lu-chandrakant/assets/Maize_China\"     # Latest dataset from Peng et al. 2023 (https://doi.org/10.1038/s41597-023-02573-6)\n",
    "    Rubber_dataset_version = 'users/wangyxtina/MapRubberPaper/rForeRub202122_perc1585DifESAdist5pxPFfinal'     # Latest dataset from Wang et al. 2023 (https://doi.org/10.1038/s41586-023-06642-z)\n",
    "\n",
    "    Cropland_dataset_version = ['users/potapovpeter/Global_cropland_2003',         # Latest dataset from Potapov et al. 2022 (https://doi.org/10.1038/s43016-021-00429-z)\n",
    "                                'users/potapovpeter/Global_cropland_2007',         # can be downloaded from https://glad.umd.edu/dataset/croplands\n",
    "                                'users/potapovpeter/Global_cropland_2011',\n",
    "                                'users/potapovpeter/Global_cropland_2015',\n",
    "                                'users/potapovpeter/Global_cropland_2019'\n",
    "                               ]\n",
    "        \n",
    "    Forest_fire_dataset_version = 'users/sashatyu/2001-2022_fire_forest_loss'          # Update Forest fire version https://glad.umd.edu/dataset/Fire_GFL\n",
    "    Dominant_forest_loss_driver_version = 'projects/lu-chandrakant/assets/Dominant_driver/Dominant_driver_2001_2022'            # Download latest Curtis et al. 2018 (https://www.science.org/doi/10.1126/science.aau3445) dataset from \n",
    "                                                                                                                                # https://data.globalforestwatch.org/ and upload to GEE\n",
    "    Forest_management_dataset_version = \"projects/lu-chandrakant/assets/Forest_Management/FML_v3-2_with-colorbar\"               # Download latest Lesiv et al. 2022 (https://www.nature.com/articles/s41597-022-01332-3) from https://zenodo.org/record/5848610 and upload to GEE\n",
    "    \n",
    "    AGB_dataset_version = ['projects/lu2-chandrakant/assets/Carbon',                   # Download the year 2000 Above ground carbon map from Harris et al. 2021 (https://www.nature.com/articles/s41558-020-00976-6) from \n",
    "                           'projects/lu3-chandrakant/assets/Carbon'                    # https://data.globalforestwatch.org/maps/e4bdbe8d6d8d4e32ace7d36a4aec7b93 and upload to GEE. \n",
    "                          ]                                                            # Since the AGB dataset is quite large, uploaded using 'geeup' package https://github.com/samapriya/geeup\n",
    "    \n",
    "    SOC_stock_dataset_version = \"projects/soilgrids-isric/ocs_mean\"                    # SoilGrid250m 2.0 (https://soilgrids.org/)\n",
    "    SOC_content_dataset_version = \"projects/soilgrids-isric/soc_mean\"                  # First 0-30m Soil Organic Carbon Stock (t/ha)\n",
    "    Bulk_density_dataset_version = \"projects/soilgrids-isric/bdod_mean\"                # Latest dataset from Poggio et al. 2021 (https://doi.org/10.5194/soil-7-217-2021) \n",
    "    Course_fragment_dataset_version = \"projects/soilgrids-isric/cfvo_mean\"\n",
    "    \n",
    "    Ecoregion_dataset_version = \"RESOLVE/ECOREGIONS/2017\"                              # To be used for BGB calculation\n",
    "    Elevation_dataset_version = 'CGIAR/SRTM90_V4'                                      # Check these datasets on Google Earth Engine\n",
    "    Precipitation_dataset_version = 'UCSB-CHG/CHIRPS/PENTAD'                           # Check these datasets on Google Earth Engine\n",
    "    Peatland_dataset_version = \"projects/lu-chandrakant/assets/Global_Peatlands\"       # Peatland dataset can be downloaded from GFW (https://data.globalforestwatch.org/datasets/gfw::global-peatlands/about)\n",
    "    \n",
    "    \n",
    "    ###################\n",
    "    # 2.1 Check value that are used for commodity classification\n",
    "    ###################\n",
    "    \n",
    "    Sugarcane_MapBiomas_reclass_value = 3221    # Code for reclassified Sugarcane from MapBiomas\n",
    "    Sugarcane_Brazil_reclass_value = 3222       # Same, but from Zheng et al.\n",
    "    \n",
    "    Soy_MapBiomas_reclass_value = 3241          # Code for reclassified Soya beans from MapBiomas\n",
    "    Soy_SouthAmerica_reclass_value = 3242       # Same, but from Song et al.\n",
    "    \n",
    "    Oilpalm_MapBiomas_reclass_value = 6121      # Code for reclassified Oil palm fruit from MapBiomas\n",
    "    Oilpalm_Plantation_reclass_value = 5121     # Same, but from Du et al.\n",
    "    Oilpalm_Global_reclass_value = 6122         # Same, but from Descals et al.\n",
    "    Oilpalm_Indonesia_reclass_value = 6123      # Same, but from Gaveau et al.\n",
    "    Oilpalm_Malaysia_reclass_value = 6124       # Same, but from Xu et al.\n",
    "    \n",
    "    Plantation_base_class_value = 5000          # Reclassified broad plantation class, added with Du et al. dataset\n",
    "    \n",
    "    Cocoa_Plantation_reclass_value = 5024       # Code for reclassified Cocoa from Du et al.\n",
    "    Cocoa_CDGandGhana_reclass_value = 6031      # Same, but from Kalischek et al.\n",
    "    \n",
    "    Coconut_Plantation_reclass_value = 5034     # Code for reclassified Coconut from Du et al.\n",
    "    Coconut_Global_reclass_value = 6041         # Same, but from Descals et al.\n",
    "    \n",
    "    Rapeseed_EuropeandCanada_reclass_value = 3301  # Code for reclassified Coconut from Han et al.\n",
    "\n",
    "    Rice_MapBiomas_reclass_value = 3261         # Code for reclassified Sugarcane from MapBiomas\n",
    "    Rice_Asia_reclass_value = 3262              # Same, but from Han et al.\n",
    "    \n",
    "    Cropland_Global_reclass_value = 3201        # Code for reclassified Cropland from Potopov et al.\n",
    "\n",
    "    Maize_China_reclass_value = 3321\n",
    "\n",
    "    Rubber_reclass_values = 6151\n",
    "\n",
    "    ###################\n",
    "    # 2.2 Check data for different coutries \n",
    "    ###################\n",
    "    Plantation_dataset_version_country = ['Argentina', 'Australia', 'Brazil', 'Cambodia', 'Cameroon', 'Chile', 'China', 'Colombia', 'Costa Rica', 'Democratic Republic of the Congo', \n",
    "                          'Ecuador', 'Gabon', 'Ghana', 'Guatemala', 'Honduras', 'India', 'Indonesia', 'Côte d\\'Ivoire', 'Japan', 'Kenya', 'Liberia', 'Malawi', 'Malaysia', 'México', \n",
    "                          'Myanmar', 'Nepal', 'New Zealand', 'Nicaragua', 'Nigeria', 'Pakistan', 'Panama', 'Papua New Guinea', 'Peru', 'Philippines', 'Rwanda', 'Solomon Islands', \n",
    "                          'South Africa', 'South Korea', 'Sri Lanka', 'Thailand', 'Uruguay', 'United States', 'Venezuela', 'Vietnam',\n",
    "                          'Åland', 'Albania', 'Andorra', 'Austria', 'Belarus', 'Belgium', 'Bosnia and Herzegovina', 'Bulgaria', 'Croatia', 'Czechia', 'Denmark', 'Estonia', \n",
    "                          'Faroe Islands', 'Finland', 'France', 'Germany', 'Greece', 'Guernsey', 'Hungary', 'Iceland', 'Ireland', 'Isle of Man', 'Italy', 'Jersey', 'Kosovo', \n",
    "                          'Latvia', 'Liechtenstein', 'Lithuania', 'Luxembourg', 'Malta', 'Moldova', 'Monaco', 'Montenegro', 'Netherlands', 'North Macedonia', 'Norway', \n",
    "                          'Poland', 'Portugal', 'Romania', 'San Marino', 'Serbia', 'Slovakia', 'Slovenia', 'Spain', 'Svalbard and Jan Mayen', 'Sweden', 'Switzerland', \n",
    "                          'Ukraine', 'United Kingdom', 'Vatican City']\n",
    "\n",
    "    MapBiomas_dataset_version_country = ['Argentina','Bolivia','Brazil','Chile','Colombia','Ecuador','French Guiana','Guyana','Paraguay','Peru','Suriname','Uruguay','Venezuela','Indonesia']\n",
    "\n",
    "    Rapeseed_dataset_version_country = ['Albania','Austria','Bulgaria','Denmark','Belarus','Estonia','Faroe Islands','Finland','France','Germany','Bosnia and Herzegovina',\n",
    "                          'Greece','Hungary','Croatia','Iceland','Ireland','Italy','Latvia','Lithuania','Malta','Moldova','Netherlands','North Macedonia','Norway',\n",
    "                          'Czechia','Poland','Portugal','Romania','Slovenia','Slovakia','Spain','Sweden','Switzerland','United Kingdom','Ukraine','Belgium',\n",
    "                          'Luxembourg','Serbia and Montenegro','Serbia','Montenegro','Åland','Andorra','Guernsey','Isle of Man','Jersey','Kosovo','Liechtenstein',\n",
    "                          'Monaco','San Marino','Svalbard and Jan Mayen','Vatican City','Turkey','United States','Canada','Chile']\n",
    "\n",
    "    Rice_dataset_version_country = ['Brunei','Myanmar','Indonesia','Cambodia','Laos','Malaysia','Philippines','Timor-Leste','Singapore','Thailand','Vietnam','China',\n",
    "                          'Japan','North Korea','South Korea','Taiwan','India']\n",
    "\n",
    "    Global_oilpalm_and_Coconut_not_in = ['Canada', 'Russia', 'United States', 'Romania', 'Croatia', 'Japan'] \n",
    "\n",
    "    \n",
    "    ####################\n",
    "    # 3. Load and pre-process datasets\n",
    "    ####################\n",
    "    \n",
    "    geometry = ee.FeatureCollection(GADM_admin_boundary_version) \n",
    "    ## Filtering the geometry to a specific country/boundary\n",
    "    if Admin_boundary in ['United States', 'Russia', 'Canada', 'Brazil']:\n",
    "        uniqueValues_ADM2CODE = geometry.filter(ee.Filter.eq('COUNTRY', Admin_boundary)).distinct('GID_2').aggregate_array('GID_2');\n",
    "        uniqueValues_ADM2CODE = uniqueValues_ADM2CODE.sort()\n",
    "        uniqueValues_ADM2CODE = uniqueValues_ADM2CODE.getInfo()\n",
    "        # Load the feature collection needed to simulate aggregation\n",
    "        geometry = geometry.filter(ee.Filter.inList('GID_2', uniqueValues_ADM2CODE)).filter(ee.Filter.eq('COUNTRY', Admin_boundary))\n",
    "        size = geometry.size().getInfo()\n",
    "        sub_folder_code = '_'\n",
    "    else:\n",
    "        uniqueValues_ADM2CODE = geometry.filter(ee.Filter.eq('COUNTRY', Admin_boundary)).distinct('COUNTRY').aggregate_array('COUNTRY');\n",
    "        # Load the feature collection needed to simulate aggregation\n",
    "        geometry = geometry.filter(ee.Filter.inList('COUNTRY', uniqueValues_ADM2CODE)).filter(ee.Filter.eq('COUNTRY', Admin_boundary))\n",
    "        size = geometry.size().getInfo()\n",
    "        sub_folder_code = ''\n",
    "    #$$size = geometry.size().getInfo()\n",
    "    print('Boundary', Admin_boundary, ':', size)\n",
    "        \n",
    "    \n",
    "    \n",
    "    ## Load tree cover and tree cover loss dataset (Hansen et al. 2013)\n",
    "    # Load the latest version of the dataset and select the variables needed to be processed\n",
    "    Hansen_data = ee.Image(Hansen_tree_cover_dataset_version).select(['treecover2000', 'loss', 'lossyear'])\n",
    "    Hansen_data = Hansen_data.clip(geometry)\n",
    "\n",
    "\n",
    "    if Admin_boundary in Plantation_dataset_version_country:\n",
    "        ## Load global plantation dataset (Du et al. 2022)\n",
    "        # Define list of image names as uplaoded in GEE\n",
    "        imageNames = [\"pYear_1\", \"pYear_2\", \"pYear_3\", \"pYear_4\", \"pYear_5\", \"pYear_6\", \"pYear_7\", \"pYear_8\", \"pYear_9\", \"pYear_10\", \"pYear_11\", \"pYear_12\", \"pYear_13\", \"pYear_14\"]\n",
    "        # Map over image names to create a list of images to create an ImageCollection for band 'b1' (i.e, plantyear)\n",
    "        images = ee.ImageCollection.fromImages(list(map(lambda name: ee.Image(Plantation_dataset_version + name)\\\n",
    "                                                        .select(['b1']).rename(name), imageNames)))\n",
    "        # Combine all images into one and concatenate all plantation bands into a single band\n",
    "        image = images.toBands().rename(imageNames)\n",
    "        concatenatedImage = image.select(imageNames)\n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        combinedImage = concatenatedImage.reduce(ee.Reducer.max())\n",
    "        # Mask regions below zero (No-data)\n",
    "        plantyear = combinedImage.updateMask(combinedImage.gt(0)).select(['max']).rename('plantyear')\n",
    "    \n",
    "        # Map over image names to create a list of images to create an ImageCollection for band 'b2' (i.e., startyear) \n",
    "        images = ee.ImageCollection.fromImages([ee.Image(Plantation_dataset_version + name)\\\n",
    "                                                .select(['b2']).rename(name) for name in imageNames])\n",
    "        image = images.toBands().rename(imageNames)\n",
    "        concatenatedImage = image.select(imageNames)\n",
    "        combinedImage = concatenatedImage.reduce(ee.Reducer.max())\n",
    "        startyear = combinedImage.updateMask(combinedImage.gt(0)).select(['max']).rename('startyear')\n",
    "    \n",
    "        #Map over image names to create a list of images to create an ImageCollection for band 'b3' (i.e., species)\n",
    "        images = ee.ImageCollection.fromImages([ee.Image(Plantation_dataset_version + name)\\\n",
    "                                                .select(['b3']).rename(name) for name in imageNames])\n",
    "        image = images.toBands().rename(imageNames)\n",
    "        concatenatedImage = image.select(imageNames)\n",
    "        combinedImage = concatenatedImage.reduce(ee.Reducer.max())\n",
    "        species = combinedImage.updateMask(combinedImage.gt(0)).select(['max']).rename('species')\n",
    "\n",
    "    \n",
    "\n",
    "    if Admin_boundary in MapBiomas_dataset_version_country:\n",
    "        # Function to retrieve datasets based on the country\n",
    "        def get_datasets_by_country(Admin_boundary):\n",
    "            datasets = {\n",
    "                'Argentina': ['Atlantic forest', 'Argentina', 'Chaco', 'Pampa'],\n",
    "                'Bolivia': ['Amazon', 'Chaco', 'Bolivia'],\n",
    "                'Brazil': ['Amazon', 'Atlantic forest', 'Brazil', 'Pampa'],\n",
    "                'Chile': ['Chile'],\n",
    "                'Colombia': ['Amazon', 'Colombia'],\n",
    "                'Ecuador': ['Amazon', 'Ecuador'],\n",
    "                'French Guiana': ['Amazon'],\n",
    "                'Guyana': ['Amazon'],\n",
    "                'Paraguay': ['Atlantic forest', 'Chaco', 'Paraguay'],\n",
    "                'Peru': ['Amazon', 'Peru'],\n",
    "                'Suriname': ['Amazon'],\n",
    "                'Uruguay': ['Uruguay', 'Pampa'],\n",
    "                'Venezuela': ['Venezuela'],\n",
    "                'Indonesia': ['Indonesia']\n",
    "            }\n",
    "            selected_keys = datasets.get(Admin_boundary, [])\n",
    "            return [MapBiomas_dataset_version[key] for key in selected_keys if key in MapBiomas_dataset_version]\n",
    "            \n",
    "        ## Load MapBiomas dataset\n",
    "        def combined_img_func(Year):\n",
    "            # First we reclassify the Mapbiomas defined classes with our classficiation. 'in_class' is Mapbiomas classification and 'reclass' is the our classification (See Supplementary Data)\n",
    "            in_class = [0,1,2,3,4,5,6,9,10,11,12,13,14,15,18,19,20,21,22,23,24,25,26,27,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,57,58,61,62,65,66]\n",
    "            reclass = [1,1000,1100,1300,2100,2100,2100,5000,2000,2100,2100,2100,3050,4000,3150,3200,3221,3100,2100,2100,600,2100,100,1,2100,700,100,2100,100,100,6121,3800,\n",
    "                                100,100,3241,3261,3200,2100,2100,2100,2100,6021,6001,3800,2100,2100,3200,3200,2100,3281,3802,2100]\n",
    "                    \n",
    "            # We analyse seperately for Indonesia and Latin America due to time-series gap. Latin America datasets are available till 2021 (when simulating this code), \n",
    "            # (contd.) while Indonesia is avaible only till 2019.   \n",
    "            imageNames = get_datasets_by_country(Admin_boundary)\n",
    "            if Admin_boundary == 'Bolivia':\n",
    "                # Map over image names to create a list of images\n",
    "                def func_pxw(imageNames):\n",
    "                    # We gapfill for year 2020 and 2021 since the Indonesia data in unavaiable then by extending 2019 dataset to end year\n",
    "                    # Change '2020' if the latest dataset becomes avaible at a later date. \n",
    "                    if (Year >= 2022) & (Year <= end_year):\n",
    "                        image = ee.Image(imageNames).select(['classification_2021'])\n",
    "                    else:\n",
    "                        image = ee.Image(imageNames).select(['classification_' + str(Year)])\n",
    "                    return image\n",
    "                images = ee.ImageCollection.fromImages([func_pxw(imageNames).rename('Class').remap(in_class, reclass, 1) for name in imageNames])\n",
    "            else:\n",
    "                # Map over image names to create a list of images\n",
    "                images = ee.ImageCollection.fromImages([ee.Image(name).select(['classification_' + str(Year)]).rename('Class').remap(in_class, reclass, 1) for name in imageNames])\n",
    "                        \n",
    "                        \n",
    "            # Combine all images into one and concatenate all plantation bands into a single band\n",
    "            image = images.toBands().rename(imageNames)\n",
    "            concatenatedImage = image.select(imageNames)\n",
    "            # Use reduce() to apply the reducer function to all bands\n",
    "            combinedImage = concatenatedImage.reduce(ee.Reducer.max());\n",
    "            #combinedImage = combinedImage.clip(geometry)#$$\n",
    "            # We take maximum so that if a pixel has multiple values\n",
    "            MapBiomas = combinedImage.select(['max']).rename('classification');\n",
    "            return MapBiomas\n",
    "        \n",
    "        \n",
    "        ## Load Soybean dataset\n",
    "        # Load the existing ImageCollection\n",
    "        soy_collection = ee.ImageCollection(Soy_dataset_version)\n",
    "    \n",
    "        def Soy_combined_img_func(Year):\n",
    "            year = Year\n",
    "            year_string = ee.String(str(year)) # convert year to string\n",
    "            # filter the soy_collection by the specified year\n",
    "            filtered_collection = soy_collection.filterDate(\n",
    "                year_string.cat('-01-01'), # start of year\n",
    "                year_string.cat('-12-31')  # end of year\n",
    "            )\n",
    "            # Use reduce() to apply the reducer function to all bands\n",
    "            combinedImage = filtered_collection.reduce(ee.Reducer.max())\n",
    "            #combinedImage = combinedImage.clip(geometry)#$$\n",
    "            combinedImage = combinedImage.updateMask(combinedImage.eq(1))\n",
    "            # Mask regions below zero (No-data)\n",
    "            soy = combinedImage.select(['b1_max']).rename('classification').multiply(Soy_SouthAmerica_reclass_value)\n",
    "            return soy\n",
    "\n",
    "\n",
    "    if Admin_boundary == 'China':\n",
    "        ## Load Maize-China (Peng et al. 2023)\n",
    "        # Specify the asset folder path where Maize features are stored\n",
    "        maize_china_list = ee.data.getList({'id': Maize_China_dataset_version})\n",
    "        maize_china_list = list(map(lambda image: image['id'].split(\"/\")[4], maize_china_list))\n",
    "        \n",
    "        # Map over image names to create an image collection for lu2-chandrakant\n",
    "        maize_china = ee.ImageCollection(list(map(lambda name: ee.Image(Maize_China_dataset_version+\"/\" + name).select(['b1']).rename('Class'), maize_china_list)))\n",
    "        \n",
    "        def maize_china_func(Year):\n",
    "            # Print the updated ImageCollection\n",
    "            if Year > 2020:\n",
    "                year = 2020\n",
    "            else:\n",
    "                year = Year\n",
    "            year_string = ee.String(str(year))  # Convert year to string\n",
    "        \n",
    "            # Filter the maize_china collection by the specified year\n",
    "            filtered_collection = maize_china.filterDate(\n",
    "                year_string.cat('-01-01T01:01:01'),  # Start of year\n",
    "                year_string.cat('-12-31T12:01:01')  # End of year\n",
    "            )\n",
    "            # Use reduce() to apply the reducer function to all bands\n",
    "            combined_image = filtered_collection.reduce(ee.Reducer.max())\n",
    "            #combined_image = combined_image.clip(geometry)#$$\n",
    "            combined_image = combined_image.updateMask(combined_image.eq(1))\n",
    "            # Mask regions below zero (No-data)\n",
    "            maize = combined_image.select(['Class_max']).rename('classification')\n",
    "            return maize.multiply(Maize_China_reclass_value)\n",
    "\n",
    "\n",
    "    if Admin_boundary in ['Côte d\\'Ivoire', 'Ghana']: \n",
    "        ## Load Cocoa dataset \n",
    "        Cocoa_map_CDG = ee.Image(Cocoa_dataset_version)\n",
    "        #Cocoa_map_CDG = Cocoa_map_CDG.clip(geometry)#$$\n",
    "        Cocoa_map_CDG_projection =  ee.Image(Cocoa_map_CDG).projection()\n",
    "\n",
    "\n",
    "    if Admin_boundary == 'Indonesia':\n",
    "        ## Load Oilpalm Indonesia dataset\n",
    "        assetList = ee.data.getList({'id': Oilpalm_Indonesia_version})\n",
    "        assetIds = [asset['id'] for asset in assetList]\n",
    "        \n",
    "        OilPalm_Indonesia_org = ee.FeatureCollection([])\n",
    "        for assetId in assetIds:\n",
    "            assetCollection = ee.FeatureCollection(assetId)\n",
    "            OilPalm_Indonesia_org = OilPalm_Indonesia_org.merge(assetCollection)\n",
    "\n",
    "\n",
    "\n",
    "    if Admin_boundary == 'Malaysia':\n",
    "        ## Load (Malaysia) Oilpalm dataset\n",
    "        # Specify the asset folder path where Oil palm features are stored\n",
    "        OilPalm_Malaysia_list = ee.data.getList({'id': Oilpalm_Malaysia_version})\n",
    "        OilPalm_Malaysia_list = [image['id'].split('/')[4] for image in OilPalm_Malaysia_list]\n",
    "        \n",
    "        # Map over image names to create an image collection for lu2-chandrakant\n",
    "        OilPalm_Malaysia = ee.ImageCollection([ee.Image(Oilpalm_Malaysia_version+ \"/\" + name).select(['b1']).rename('Class') for name in OilPalm_Malaysia_list])\n",
    "        \n",
    "        \n",
    "        def OilPalm_Malaysia_func(Year):\n",
    "            # Print the updated ImageCollection\n",
    "            if Year > 2018:\n",
    "                year = 2018\n",
    "            else:\n",
    "                year = Year\n",
    "        \n",
    "            yearString = ee.String(str(year))  # convert year to string\n",
    "            # filter the OilPalm_Malaysia collection by the specified year\n",
    "            filteredCollection = OilPalm_Malaysia.filterDate(\n",
    "                yearString.cat('-01-01'),  # start of year\n",
    "                yearString.cat('-12-31')  # end of year\n",
    "            )\n",
    "        \n",
    "            # Use reduce() to apply the reducer function to all bands\n",
    "            combinedImage = filteredCollection.reduce(ee.Reducer.max())\n",
    "            #combinedImage = combinedImage.clip(geometry)#$$\n",
    "            combinedImage = combinedImage.updateMask(combinedImage.eq(1))\n",
    "        \n",
    "            # Mask regions below zero (No-data)\n",
    "            OilPalm = combinedImage.select(['Class_max']).rename('classification')\n",
    "            return OilPalm.multiply(Oilpalm_Malaysia_reclass_value)\n",
    "\n",
    "\n",
    "    if Admin_boundary not in Global_oilpalm_and_Coconut_not_in:\n",
    "        ## Load (Global) Oilpalm dataset\n",
    "        OilPalm_Global = ee.ImageCollection(Oilpalm_Global_version)\n",
    "        OilPalm_Global_projection =  ee.Image(OilPalm_Global.first()).projection()      # Will be needed later for ee.reduceResolution\n",
    "        OilPalm_Global = OilPalm_Global.reduce(ee.Reducer.min())\n",
    "        #OilPalm_Global = OilPalm_Global.clip(geometry)#$$\n",
    "        OilPalm_Global = OilPalm_Global.updateMask(OilPalm_Global.neq(3))\n",
    "        OilPalm_Global = OilPalm_Global.select(['classification_min']).rename(['Class'])\n",
    "    \n",
    "    \n",
    "    \n",
    "        ## Load Coconut dataset\n",
    "        Coconut_list = ee.data.getList({'id': Coconut_dataset_version})\n",
    "        Coconut_list = [((image['id']).split(\"/\")[4]) for image in Coconut_list]\n",
    "        Coconut = ee.ImageCollection.fromImages([ee.Image(Coconut_dataset_version+\"/\" + name).select(['b1']).rename('Class') for name in Coconut_list])\n",
    "        Coconut_projection =  ee.Image(Coconut.first()).projection()\n",
    "        Coconut = Coconut.reduce(ee.Reducer.min()).rename('Class')\n",
    "        #Coconut = Coconut.clip(geometry)#$$\n",
    "        Coconut = Coconut.updateMask(Coconut.eq(1))\n",
    "    \n",
    "\n",
    "    \"\"\"\n",
    "    ## Load Rubber dataset (Wang et al. 2023)\n",
    "    Rubber = ee.Image(Rubber_dataset_version)\n",
    "    #Rubber = Rubber.clip(geometry)#$$\n",
    "    Rubber_projection =  ee.Image(Rubber).projection()\n",
    "    Rubber = Rubber.updateMask(Rubber.eq(2)) # 1-Forest and 2-Rubber\n",
    "    Rubber = Rubber.where(Rubber, 1)\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    ## Load Cropland datasets\n",
    "    Cropland_2000_03 = ee.ImageCollection(Cropland_dataset_version[0]).reduce(ee.Reducer.max()).rename('Class')    \n",
    "    Cropland_2004_07 = ee.ImageCollection(Cropland_dataset_version[1]).reduce(ee.Reducer.max()).rename('Class')    \n",
    "    Cropland_2008_11 = ee.ImageCollection(Cropland_dataset_version[2]).reduce(ee.Reducer.max()).rename('Class')    \n",
    "    Cropland_2012_15 = ee.ImageCollection(Cropland_dataset_version[3]).reduce(ee.Reducer.max()).rename('Class')    \n",
    "    Cropland_2016_19 = ee.ImageCollection(Cropland_dataset_version[4]).reduce(ee.Reducer.max()).rename('Class')\n",
    "\n",
    "\n",
    "    if Admin_boundary == 'Brazil':\n",
    "        ## Load Sugarcane Brazil dataset\n",
    "        Sugarcane_list = ee.data.getList({'id': Sugarcane_dataset_version})\n",
    "        Sugarcane_list = [((image['id']).split(\"/\")[4]) for image in Sugarcane_list]\n",
    "        Sugarcane_Brazil = ee.ImageCollection.fromImages([ee.Image(Sugarcane_dataset_version+\"/\" + name).select(['b1']).rename('Class') for name in Sugarcane_list])\n",
    "        Sugarcane_Brazil = Sugarcane_Brazil.reduce(ee.Reducer.max()).rename('Class')\n",
    "        #Sugarcane_Brazil = Sugarcane_Brazil.clip(geometry)#$$\n",
    "        Sugarcane_Brazil = Sugarcane_Brazil.updateMask(Sugarcane_Brazil.eq(1))\n",
    "\n",
    "    \n",
    "\n",
    "    if Admin_boundary in Rapeseed_dataset_version_country:\n",
    "        ## Load Rapeseed dataset\n",
    "        Rapeseed_list = ee.data.getList({'id': Rapeseed_dataset_version})\n",
    "        Rapeseed_list = [((image['id']).split(\"/\")[4]) for image in Rapeseed_list]\n",
    "        Rapeseed = ee.ImageCollection.fromImages([ee.Image(Rapeseed_dataset_version+\"/\" + name).select(['b1']).rename('Class') for name in Rapeseed_list])\n",
    "        Rapeseed_projection =  ee.Image(Rapeseed.first()).projection()\n",
    "        Rapeseed = Rapeseed.reduce(ee.Reducer.max()).rename('Class')\n",
    "        #Rapeseed = Rapeseed.clip(geometry)#$$\n",
    "        Rapeseed = Rapeseed.updateMask(Rapeseed.eq(1))\n",
    "\n",
    "\n",
    "    if Admin_boundary in Rice_dataset_version_country:\n",
    "        ## Load Paddy rice dataset\n",
    "        Paddy_rice_list = ee.data.getList({'id': Rice_dataset_version})\n",
    "        Paddy_rice_list = [((image['id']).split(\"/\")[4]) for image in Paddy_rice_list]\n",
    "        Paddy_rice = ee.ImageCollection.fromImages([ee.Image(Rice_dataset_version+\"/\" + name).select(['b1']).rename('Class') for name in Paddy_rice_list])\n",
    "        Paddy_rice_projection =  ee.Image(Paddy_rice.first()).projection()\n",
    "        Paddy_rice = Paddy_rice.reduce(ee.Reducer.max()).rename('Class')\n",
    "        #Paddy_rice = Paddy_rice.clip(geometry)#$$\n",
    "        Paddy_rice = Paddy_rice.updateMask(Paddy_rice.eq(1))\n",
    "\n",
    "\n",
    "\n",
    "    ## Load Forest management data\n",
    "    ForestManagement = ee.Image(Forest_management_dataset_version);\n",
    "    #ForestManagement = ForestManagement.clip(geometry)#$$\n",
    "    \n",
    "    \n",
    "    ## Load forest loss due to fire\n",
    "    Forest_fire = ee.ImageCollection(Forest_fire_dataset_version)\n",
    "    # 1-No fire, 2- Low certainity, 3-Medium certainity, 4-High certainity and 5-'all forest loss due to fire in Africa'\n",
    "    in_class = [1, 2, 3, 4, 5]\n",
    "    reclass = [1, 1, 250, 250, 250]\n",
    "    Forest_fire = Forest_fire.reduce(ee.Reducer.min()).select(['b1_min']).rename('classification')\n",
    "    #Forest_fire = Forest_fire.clip(geometry)#$$\n",
    "    Forest_fire = Forest_fire.remap(in_class, reclass, 1)\n",
    "\n",
    "    \n",
    "    \n",
    "    ## Load Dominant tree cover loss driver\n",
    "    dominant_driver = ee.Image(Dominant_forest_loss_driver_version)\n",
    "    #dominant_driver = dominant_driver.clip(geometry)#$$\n",
    "    in_class = [1, 2, 3, 4, 5] \n",
    "    # These classes mean 1-Commodity-driven deforestation, 2-Shifting agriculture, 3-Forestry, 4-Wildfire, and 5-Urbanization\n",
    "    # Reclassify them to our classification (See Supplementary Data)\n",
    "    reclass = [3000, 3000, 500, 200, 600] \n",
    "    dominant_driver = dominant_driver.remap(in_class, reclass, 1)\n",
    "\n",
    "    \n",
    "    \n",
    "    ## Load Above ground biomass dataset (Harris et al. 2018)\n",
    "    if simulation in ['All', 'AGB', 'BGB']:\n",
    "        # Get the list of images in the lu2-chandrakant asset (i.e., uploaded to GEE)\n",
    "        lu2_asset_list = ee.data.getList({'id': AGB_dataset_version[0]})\n",
    "        lu2_image_names = [image['id'].split('/')[-1] for image in lu2_asset_list]\n",
    "        # Map over image names to create an image collection for lu2-chandrakant\n",
    "        lu2_images = ee.ImageCollection.fromImages([ee.Image(AGB_dataset_version[0]+f'/{name}').select(['b1']).rename(name) for name in lu2_image_names])\n",
    "        # Get the list of images in the lu3-chandrakant asset\n",
    "        lu3_asset_list = ee.data.getList({'id': AGB_dataset_version[1]})\n",
    "        lu3_image_names = [image['id'].split('/')[-1] for image in lu3_asset_list]\n",
    "        # Map over image names to create an image collection for lu3-chandrakant\n",
    "        lu3_images = ee.ImageCollection.fromImages([ee.Image(AGB_dataset_version[1]+f'/{name}').select(['b1']).rename(name) for name in lu3_image_names])\n",
    "    \n",
    "        # Merge the two image collections into one\n",
    "        all_images = lu2_images.merge(lu3_images)\n",
    "        # Combine all images into one and concatenate all plantation bands into a single band\n",
    "        image = all_images.toBands().rename(lu2_image_names + lu3_image_names)\n",
    "        concatenated_image = image.select(lu2_image_names + lu3_image_names)\n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        combined_image = concatenated_image.reduce(ee.Reducer.max())\n",
    "        #combined_image = combined_image.clip(geometry)#$$\n",
    "        # Mask regions below zero (No-data)\n",
    "        AGB = combined_image.updateMask(combined_image.gte(0)).select(['max']).rename('MG_px-1')\n",
    "\n",
    "\n",
    "        ## Load elevation data and apply mask\n",
    "        elevation = ee.Image(Elevation_dataset_version).select('elevation')\n",
    "        #elevation = elevation.clip(geometry)#$$\n",
    "        ## Load and filter precipitation data\n",
    "        precipitation = ee.ImageCollection(Precipitation_dataset_version) \\\n",
    "            .filter(ee.Filter.date('2000-01-01', '2000-12-31')) \\\n",
    "            .select('precipitation') \\\n",
    "            .reduce(ee.Reducer.sum())\n",
    "        #precipitation = precipitation.clip(geometry)#$$\n",
    "    \n",
    "\n",
    "    \n",
    "    ## Load SOC database\n",
    "    if simulation in ['All', 'SOC 0-30', 'SOC 30-100']:\n",
    "        # Function to calculate SOC for a given depth range\n",
    "        # Combine SOC layers for different depth ranges\n",
    "        SOC_0_30 = ee.Image(SOC_stock_dataset_version).select(\"ocs_0-30cm_mean\");    # SOC (t/ha) is already derived from SoilGrid\n",
    "        #SOC_0_30 = SOC_0_30.clip(geometry)#$$\n",
    "        def calculate_SOC(depthRange1, depthRange2):\n",
    "            soc = ee.Image(SOC_content_dataset_version).select(\"soc_\" + str(depthRange1) + '-' + str(depthRange2) + \"cm_mean\")\n",
    "            bd = ee.Image(Bulk_density_dataset_version).select(\"bdod_\" + str(depthRange1) + '-' + str(depthRange2) + \"cm_mean\")\n",
    "            cf = ee.Image(Course_fragment_dataset_version).select(\"cfvo_\" + str(depthRange1) + '-' + str(depthRange2) + \"cm_mean\")\n",
    "            cf = (cf.multiply(1e-3).subtract(1)).multiply(-1)\n",
    "            soc_t_per_ha = (soc.multiply(bd).multiply(1e-5).multiply(depthRange2-depthRange1)).multiply(cf).rename('SOC t ha-1');\n",
    "            #soc_t_per_ha = soc_t_per_ha.clip(geometry)#$$\n",
    "            return soc_t_per_ha\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    ## Load Global Peatland dataset\n",
    "    if simulation in ['All', 'PEATLAND', 'SOC 0-30', 'SOC 30-100']:\n",
    "        # Get a list of images in the asset collection\n",
    "        Peatland_list = ee.data.getList({'id': Peatland_dataset_version})\n",
    "        # Extract image names from the list\n",
    "        Peatland_list = [image['id'].split('/')[-1] for image in Peatland_list]\n",
    "        # Map over image names to create an image collection\n",
    "        Peatland = ee.ImageCollection([ee.Image(f\"{Peatland_dataset_version}/{name}\").select(['b1']).rename('Class') for name in Peatland_list])\n",
    "        # Reduce the image collection by taking the maximum pixel value\n",
    "        Peatland = Peatland.reduce(ee.Reducer.max()).rename('Class')\n",
    "        #Peatland = Peatland.clip(geometry)#$$\n",
    "        # Select pixels with a value of 1 (representing Peatland)\n",
    "        Peatland = Peatland.updateMask(Peatland.eq(1))\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "    ####################\n",
    "    # 4. Analysing dataset\n",
    "    ####################\n",
    "    \n",
    "    ## Analysing forest loss dataset\n",
    "    # List sequence required for processing. Here '1' refers to 2001.\n",
    "    #years = ee.List.sequence(1, end_year-2000)\n",
    "    # Extract the scale of the Hansen variable. This will be needed for resampling other datasets.\n",
    "    Hansen_scale = Hansen_data.projection().nominalScale()\n",
    "    # Calculate the tree cover area\n",
    "    Hansen_treeCover = Hansen_data.select(['treecover2000'])\n",
    "    Hansen_loss = Hansen_data.select(['loss'])\n",
    "    Hansen_lossyear = Hansen_data.select(['lossyear'])\n",
    "\n",
    "    # Hansen projection\n",
    "    Hansen_projection = ee.Image(Hansen_loss).projection()\n",
    "    \n",
    "    # Restricting forest loss to Forest_threshold\n",
    "    Hansen_treeCover = Hansen_treeCover.updateMask(Hansen_treeCover.gte(Forest_threshold))\n",
    "    TC_mask = Hansen_loss.gt(0).And(Hansen_treeCover.gt(0))\n",
    "    #Updating mask\n",
    "    Hansen_data = Hansen_data.updateMask(TC_mask)\n",
    "    Hansen_loss = Hansen_loss.updateMask(TC_mask)\n",
    "\n",
    "    # Hansen's data has the treecover2000 layer ranging from 0-100.\n",
    "    Hansen_treeCover = Hansen_treeCover.divide(100)\n",
    "    # It needs to be multiplied by 10^(-4) to convert the areas from 'm2' to 'ha'.\n",
    "    #area_HansenLoss = Hansen_loss.gt(0).multiply(ee.Image.pixelArea()).divide(1e4).select([0], [\"arealoss\"])\n",
    "\n",
    "    \n",
    "    if Admin_boundary in MapBiomas_dataset_version_country:\n",
    "        ## Analysing MapBiomas dataset\n",
    "        # Define the years list\n",
    "        years = [year for year in range(2001, end_year+1)]\n",
    "        # Define a function to get MapBiomas data for a given year and apply a mask (i.e., mask removes all regions where forest loss isn't happening)\n",
    "        def getMapBiomas(year):\n",
    "            # Load the MapBiomas data for the given year\n",
    "            MapBiomas = combined_img_func(year).updateMask(TC_mask)\n",
    "            # Return the masked MapBiomas image with the year as a property\n",
    "            return MapBiomas.set('year', year)\n",
    "    \n",
    "        # Map the getMapBiomas function over the years array to create an ImageCollection\n",
    "        MapBiomas_collection = ee.ImageCollection.fromImages(\n",
    "            list(map(getMapBiomas, years))\n",
    "        )\n",
    "        \n",
    "        MapBiomasyear = ee.List.sequence(2001, end_year, 1)\n",
    "        # Sequencing over only one year (i.e., forest loss year) and extracting values to MapBiomas pixel\n",
    "        \"\"\"\n",
    "        # Map over image names to create a list of images\n",
    "        images = ee.ImageCollection.fromImages(\n",
    "            MapBiomasyear.map(\n",
    "                lambda name: ee.Image(\n",
    "                    MapBiomas_collection.filterMetadata('year', 'equals', ee.Number(name)).first())\n",
    "                .select('classification')\n",
    "                .updateMask(Hansen_lossyear.eq(ee.Number(name).subtract(2000)))\n",
    "                .rename('classification'))\n",
    "        )\n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        MapBiomas = images.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "        \"\"\"\n",
    "        \n",
    "        # Sequencing over only four year (i.e., forest loss year + three proceeding years) and extracting maximum pixel value (i.e., commodity-driven)\n",
    "        # This way, plantation will be preferred over perrinial crops, and perrinial crops preferred over temperory crops\n",
    "        # Map over image names to create a list of images\n",
    "        images = ee.ImageCollection.fromImages(MapBiomasyear.map(lambda name: \n",
    "            ee.ImageCollection(MapBiomas_collection).filter(ee.Filter.inList('year', ee.List.sequence(ee.Number(name), ee.Number(name).add(3), 1)))\\\n",
    "                .reduce(ee.Reducer.max()).select(['classification_max'])\\\n",
    "                .rename('classification')\\\n",
    "                .updateMask(Hansen_lossyear.eq(ee.Number(name).subtract(2000)))\\\n",
    "                .rename('classification')\n",
    "        ))\n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        MapBiomas = images.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "        # Matching our pixels with those observed for year 2000\n",
    "        MapBiomas_2000 = combined_img_func(2000).updateMask(TC_mask)\n",
    "        # If MapBiomas 2000 is equal to MapBiomas(year), the pixel value is multiplied by -1 to suggest pre-2000's disturbance\n",
    "        MapBiomas = MapBiomas.where(MapBiomas.eq(MapBiomas_2000), MapBiomas.multiply(-1))\n",
    "        # Adding this line to remove negative 1 values (since 1 is classified as unknown)\n",
    "        MapBiomas = MapBiomas.where(MapBiomas.eq(-1), 1)\n",
    "    \n",
    "        \n",
    "        \n",
    "        ## Analysing Soybean dataset    \n",
    "        years = [year for year in range(2001, end_year+1)]\n",
    "        # Create an empty ImageCollection to store Soya beans annual dataset\n",
    "        Soy_collection = ee.ImageCollection([])\n",
    "        # Define a function to get Soya beans data for a given year and apply a mask\n",
    "        def getSoy(year):\n",
    "            # Load the Soya beans data for the given year\n",
    "            Soy = ee.Image(Soy_combined_img_func(year)).updateMask(TC_mask)\n",
    "            # Return the masked Soya beans image with the year as a property\n",
    "            return Soy.set('year', year)\n",
    "        # Map the getSoy function over the years array to create an ImageCollection\n",
    "        Soy_collection = ee.ImageCollection.fromImages(\n",
    "            list(map(getSoy, years))\n",
    "        )\n",
    "        \n",
    "        Soyyear = ee.List.sequence(2001, end_year, 1)\n",
    "        # The code below is sequncing over only one year (i.e., forest loss year)\n",
    "        \"\"\"\n",
    "        # The year of forest loss was same as Hansen forest loss\n",
    "        def get_masked_Soy(name):\n",
    "            year = ee.Number(name) # Convert name to a number\n",
    "            image = Soy_collection.filterMetadata('year', 'equals', year).first()\n",
    "            maskedImage = ee.Image(image).select('classification').updateMask(Hansen_lossyear.eq(year.subtract(2000)))\n",
    "            return maskedImage.rename('classification')\n",
    "    \n",
    "        # Map the get_masked_Soy function over the years array to create an ImageCollection\n",
    "        images = ee.ImageCollection.fromImages(Soyyear.map(get_masked_Soy))\n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        Soy = images.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "        \"\"\"\n",
    "        \n",
    "        # The code below is sequncing over four year (i.e., forest loss year + three proceeding years)\n",
    "        # Map the Soy_collection function over the years array to create an ImageCollection\n",
    "        def get_four_years_Soy(name):\n",
    "            year = ee.Number(name) # Convert name to a number\n",
    "            years = ee.List.sequence(year, year.add(3), 1)\n",
    "            filter = ee.Filter.inList('year', years)\n",
    "            image = Soy_collection.filter(filter)\n",
    "            image = image.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "            maskedImage = ee.Image(image).select('classification').updateMask(Hansen_lossyear.eq(year.subtract(2000)))\n",
    "            return maskedImage.rename('classification')\n",
    "        images = ee.ImageCollection.fromImages(Soyyear.map(get_four_years_Soy))\n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        Soy_four = images.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "\n",
    "    \n",
    "\n",
    "    if Admin_boundary == 'China':\n",
    "        # Analyse Maize-China dataset\n",
    "        # Define a list of years\n",
    "        years = [year for year in range(2001, end_year+1)]\n",
    "        \n",
    "        # Define a function to get Maize-China data for a given year and apply a mask\n",
    "        def get_maize(year):\n",
    "            # Load the Maize-China data for the given year\n",
    "            maize = maize_china_func(year).updateMask(TC_mask)\n",
    "            # Return the masked Maize-China image with the year as a property\n",
    "            return maize.set('year', year)\n",
    "        \n",
    "        # Map the get_maize function over the years list to create an ImageCollection\n",
    "        maize_china_collection = ee.ImageCollection(list(map(get_maize, years)))\n",
    "        \n",
    "        # Define a list of years for processing\n",
    "        maize_years = ee.List.sequence(2001, end_year, 1)\n",
    "    \n",
    "        # Map the get_maize function over the maize_years list to create an ImageCollection\n",
    "        def get_four_years_Maize(name):\n",
    "            year = ee.Number(name) # Convert name to a number\n",
    "            years = ee.List.sequence(year, year.add(3), 1)\n",
    "            filter = ee.Filter.inList('year', years)\n",
    "            image = maize_china_collection.filter(filter)\n",
    "            image = image.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "            maskedImage = ee.Image(image).select('classification').updateMask(Hansen_lossyear.eq(year.subtract(2000)))\n",
    "            return maskedImage.rename('classification')\n",
    "        images = ee.ImageCollection.fromImages(maize_years.map(get_four_years_Maize))\n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        Maize_China_four = images.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "\n",
    "        \n",
    "\n",
    "    if Admin_boundary == 'Indonesia':\n",
    "        ## Analyse Oilpalm Indonesia dataset    \n",
    "        # Filter out non-forest and smallholder plantations in year 2000\n",
    "        OilPalm_Indonesia = OilPalm_Indonesia_org.filter(\n",
    "          ee.Filter.And(\n",
    "            ee.Filter.neq('F2000', 'IOPP'),\n",
    "            ee.Filter.neq('F2000', 'Smallholder')\n",
    "          ))\n",
    "        OilPalm_Indonesia_pre2000 = OilPalm_Indonesia_org.filter(\n",
    "          ee.Filter.Or(\n",
    "            ee.Filter.eq('F2000', 'IOPP'),\n",
    "            ee.Filter.eq('F2000', 'Smallholder')\n",
    "          ))\n",
    "        \n",
    "        # Add the new property to each feature in the feature collection\n",
    "        def set_reclass_code(feature):\n",
    "            return feature.set('Reclass Code', Oilpalm_Indonesia_reclass_value)\n",
    "        \n",
    "        OilPalm_Indonesia = OilPalm_Indonesia.map(set_reclass_code)\n",
    "        OilPalm_Indonesia_pre2000 = OilPalm_Indonesia_pre2000.map(set_reclass_code)\n",
    "        OilPalm_Indonesia_Image = OilPalm_Indonesia.reduceToImage(\n",
    "          properties=['Reclass Code'],        # replace 'property_name' with the name of the property in your shapefile that contains the value you want to use\n",
    "          reducer=ee.Reducer.sum()\n",
    "        )\n",
    "        OilPalm_Indonesia_pre2000_Image = OilPalm_Indonesia_pre2000.reduceToImage(\n",
    "          properties=['Reclass Code'],\n",
    "          reducer=ee.Reducer.sum()\n",
    "        )\n",
    "        # Clip and update masks\n",
    "        #OilPalm_Indonesia_Image = OilPalm_Indonesia_Image.clip(geometry)#$$\n",
    "        #OilPalm_Indonesia_pre2000_Image = OilPalm_Indonesia_pre2000_Image.clip(geometry)#$$\n",
    "    \n",
    "        OilPalm_Indonesia_Image = OilPalm_Indonesia_Image.updateMask(TC_mask)\n",
    "        OilPalm_Indonesia_pre2000_Image = OilPalm_Indonesia_pre2000_Image.updateMask(TC_mask)\n",
    "\n",
    "\n",
    "    \n",
    "    ## Analyze Cropland dataset\n",
    "    # Selecting the dataset within the appropriate year range\n",
    "    Cropland_2000_03 = Cropland_2000_03.updateMask(Cropland_2000_03.eq(1).And(Hansen_lossyear.gte(1)).And(Hansen_lossyear.lte(3)))\n",
    "    Cropland_2004_07 = Cropland_2004_07.updateMask(Cropland_2004_07.eq(1).And(Hansen_lossyear.gte(1)).And(Hansen_lossyear.lte(7)))\n",
    "    Cropland_2008_11 = Cropland_2008_11.updateMask(Cropland_2008_11.eq(1).And(Hansen_lossyear.gte(5)).And(Hansen_lossyear.lte(11)))\n",
    "    Cropland_2012_15 = Cropland_2012_15.updateMask(Cropland_2012_15.eq(1).And(Hansen_lossyear.gte(8)).And(Hansen_lossyear.lte(15)))\n",
    "    Cropland_2016_19 = Cropland_2016_19.updateMask(Cropland_2016_19.eq(1).And(Hansen_lossyear.gte(12)).And(Hansen_lossyear.lte(19)))\n",
    "    \n",
    "    # Create an ImageCollection with the cropland expansion datasets\n",
    "    combinedCropland = ee.ImageCollection([\n",
    "      Cropland_2000_03.updateMask(TC_mask),\n",
    "      Cropland_2004_07.updateMask(TC_mask),\n",
    "      Cropland_2008_11.updateMask(TC_mask),\n",
    "      Cropland_2012_15.updateMask(TC_mask),\n",
    "      Cropland_2016_19.updateMask(TC_mask)\n",
    "    ])\n",
    "    # Mosaic the ImageCollection to create a single image\n",
    "    combinedCropland = combinedCropland.mosaic()\n",
    "    # Reduce and rename the combined image\n",
    "    combinedCropland = combinedCropland.reduce(ee.Reducer.max()).rename('Class')\n",
    "    #combinedCropland = combinedCropland.clip(geometry)#$$\n",
    "    combinedCropland = combinedCropland.multiply(Cropland_Global_reclass_value)\n",
    "    \n",
    "\n",
    "\n",
    "    if Admin_boundary == 'Brazil':\n",
    "        ## Analyze Sugarcane Brazil dataset\n",
    "        Sugarcane_Brazil = Sugarcane_Brazil.updateMask(TC_mask)\n",
    "        Sugarcane_Brazil = Sugarcane_Brazil.where(Sugarcane_Brazil.gt(0), Sugarcane_Brazil_reclass_value)\n",
    "\n",
    "    \n",
    "\n",
    "    if Admin_boundary in Rapeseed_dataset_version_country:\n",
    "        ## Analyze Rapeseed dataset\n",
    "        Rapeseed = Rapeseed.updateMask(TC_mask)\n",
    "        # Set a default projection\n",
    "        Rapeseed = Rapeseed.setDefaultProjection(Rapeseed_projection)\n",
    "        Rapeseed = Rapeseed.reduceResolution(reducer = ee.Reducer.count(), bestEffort = True, maxPixels = 3*3).reproject(Hansen_projection)\n",
    "        Rapeseed = Rapeseed.updateMask(Rapeseed.gt(4))       # Reduce resolution based on majority interpolation \n",
    "        Rapeseed = Rapeseed.where(Rapeseed.gt(0), Rapeseed_EuropeandCanada_reclass_value)\n",
    "\n",
    "\n",
    "\n",
    "    if Admin_boundary in Rice_dataset_version_country:\n",
    "        ## Analyze Paddy rice dataset\n",
    "        Paddy_rice = Paddy_rice.updateMask(TC_mask)\n",
    "        Paddy_rice = Paddy_rice.setDefaultProjection(Paddy_rice_projection)\n",
    "        Paddy_rice = Paddy_rice.reduceResolution(reducer = ee.Reducer.count(), bestEffort = True, maxPixels = 3*3).reproject(Hansen_projection)\n",
    "        Paddy_rice = Paddy_rice.updateMask(Paddy_rice.gt(4))       # Reduce resolution based on majority interpolation \n",
    "        Paddy_rice = Paddy_rice.where(Paddy_rice.gt(0), Rice_Asia_reclass_value)\n",
    "    \n",
    "    \n",
    "\n",
    "    if Admin_boundary in Plantation_dataset_version_country:\n",
    "        ## Analysing Plantation dataset\n",
    "        # A value of '1981' means that the planting year was before 1982\n",
    "        plantyear = plantyear.updateMask(TC_mask)\n",
    "        startyear = startyear.updateMask(TC_mask)\n",
    "        startyear = startyear.updateMask(plantyear.gte(1))\n",
    "        startyear = startyear.updateMask(plantyear.gt(startyear)) # Plant year shouldn't be less than start year\n",
    "        species = species.updateMask(TC_mask).updateMask(startyear)\n",
    "        \n",
    "        # Seperating pre- and post-2000's deforestation using startyear\n",
    "        FL_att_to_new_plantations = startyear.updateMask(startyear.gt(Plantation_threshold_year))\n",
    "        FL_att_to_old_plantations = startyear.updateMask(startyear.lte(Plantation_threshold_year))\n",
    "        FL_att_to_new_plantations_species = species.updateMask(FL_att_to_new_plantations.gt(0))\n",
    "        FL_att_to_old_plantations_species = species.updateMask(FL_att_to_old_plantations.gt(0))\n",
    "    \n",
    "        # Maching with the classification of MapBiomas-reclass\n",
    "        FL_att_to_old_plantations_species = FL_att_to_old_plantations_species.add(Plantation_base_class_value)\n",
    "        FL_att_to_new_plantations_species = FL_att_to_new_plantations_species.add(Plantation_base_class_value)\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## Complimenting Curtis et al with Forest management dataset to screen pre- and post-2000 deforestion byu removing already managed forests\n",
    "    ForestManagement = ForestManagement.updateMask(TC_mask)\n",
    "    ManagedForests_org = ForestManagement.updateMask(ForestManagement.eq(20)                          # 20: Naturally regenerating forest with signs of management, e.g., logging, clear cuts etc;\n",
    "                                                 .Or(ForestManagement.eq(31))                         # 31: Planted forests;\n",
    "                                                 .Or(ForestManagement.eq(32))                         # 32: Plantation forests (rotation time up to 15 years);\n",
    "                                                 .Or(ForestManagement.eq(40))                         # 40: Oil palm plantations;\n",
    "                                                 .Or(ForestManagement.eq(53)))                        # 53: Agroforestry;\n",
    "\n",
    "    # Check if Admin_boundary is in the country list (i.e., Du et al. Forest plantatation data is preferred)\n",
    "    if Admin_boundary in Plantation_dataset_version_country:\n",
    "        ManagedForests = FL_att_to_old_plantations_species\n",
    "    else: #Otherwise Managed Forest dataset is preferred \n",
    "        ManagedForests = ManagedForests_org\n",
    "\n",
    "\n",
    "\n",
    "    if Admin_boundary in ['Côte d\\'Ivoire', 'Ghana']: \n",
    "        ## Analyze Cocoa dataset\n",
    "        Cocoa_map_CDG = Cocoa_map_CDG.updateMask(TC_mask)\n",
    "        Cocoa_map_CDG = Cocoa_map_CDG.setDefaultProjection(Cocoa_map_CDG_projection)\n",
    "        Cocoa_map_CDG = Cocoa_map_CDG.reduceResolution(reducer = ee.Reducer.count(), bestEffort = True, maxPixels = 3*3).reproject(Hansen_projection)\n",
    "        Cocoa_map_CDG = Cocoa_map_CDG.updateMask(Cocoa_map_CDG.gt(4))       # Reduce resolution based on majority interpolation \n",
    "        Cocoa_map_CDG = Cocoa_map_CDG.where(Cocoa_map_CDG.gt(0), Cocoa_CDGandGhana_reclass_value)\n",
    "        Cocoa_map_CDG = Cocoa_map_CDG.where(ManagedForests.gt(0).And(Cocoa_map_CDG.gt(1)), Cocoa_map_CDG.multiply(-1))\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "    if Admin_boundary == 'Malaysia': \n",
    "        ## Analyze (Malaysia) Oilpalm dataset\n",
    "        years = [year for year in range(2001, end_year+1)]\n",
    "        # Define a function to get OilPalm data for a given year and apply a mask\n",
    "        def getOilPalm(year):\n",
    "            # Load the OilPalm data for the given year\n",
    "            OilPalm = OilPalm_Malaysia_func(year).updateMask(TC_mask)\n",
    "            # Return the masked OilPalm image with the year as a property\n",
    "            return OilPalm.set('year', year)\n",
    "        \n",
    "        # Map the getOilPalm function over the years list to create an ImageCollection\n",
    "        OilPalm_Malaysia_collection = ee.ImageCollection.fromImages(\n",
    "            list(map(getOilPalm, years))\n",
    "        )\n",
    "        \n",
    "        # Create a list of years for processing\n",
    "        OilPalmyear = ee.List.sequence(2001, end_year, 1)\n",
    "        \n",
    "        # Map the getOilPalm function over the years list to create an ImageCollection\n",
    "        def get_four_years_OilPalm(name):\n",
    "            year = ee.Number(name) # Convert name to a number\n",
    "            years = ee.List.sequence(year, year.add(3), 1)\n",
    "            filter = ee.Filter.inList('year', years)\n",
    "            image = OilPalm_Malaysia_collection.filter(filter)\n",
    "            image = image.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "            maskedImage = ee.Image(image).select('classification').updateMask(Hansen_lossyear.eq(year.subtract(2000)))\n",
    "            return maskedImage.rename('classification')\n",
    "        images = ee.ImageCollection.fromImages(OilPalmyear.map(get_four_years_OilPalm))\n",
    "    \n",
    "        # Use reduce() to apply the reducer function to all bands\n",
    "        OilPalm_Malaysia_four = images.reduce(ee.Reducer.max()).select(['classification_max']).rename('classification')\n",
    "        OilPalm_Malaysia_four = OilPalm_Malaysia_four.where(ManagedForests.gt(0).And(OilPalm_Malaysia_four.gt(0)), OilPalm_Malaysia_four.multiply(-1))\n",
    "\n",
    "    \n",
    "    if Admin_boundary not in Global_oilpalm_and_Coconut_not_in:\n",
    "        ## Analyze (Global) Oilpalm dataset\n",
    "        OilPalm_Global = OilPalm_Global.updateMask(TC_mask)\n",
    "        OilPalm_Global = OilPalm_Global.setDefaultProjection(OilPalm_Global_projection)\n",
    "        OilPalm_Global = OilPalm_Global.reduceResolution(reducer = ee.Reducer.count(), bestEffort = True, maxPixels = 3*3).reproject(Hansen_projection)\n",
    "        OilPalm_Global = OilPalm_Global.updateMask(OilPalm_Global.gt(4))       # Reduce resolution based on majority interpolation \n",
    "        OilPalm_Global = OilPalm_Global.where(OilPalm_Global.gt(0), Oilpalm_Global_reclass_value)\n",
    "        OilPalm_Global = OilPalm_Global.where(ManagedForests.gt(0).And(OilPalm_Global.gt(0)), OilPalm_Global.multiply(-1))\n",
    "        \n",
    "        \n",
    "        ## Analyze Coconut dataset\n",
    "        Coconut = Coconut.updateMask(TC_mask)\n",
    "        Coconut = Coconut.setDefaultProjection(Coconut_projection)\n",
    "        Coconut = Coconut.reduceResolution(reducer = ee.Reducer.count(), bestEffort = True, maxPixels = 3).reproject(Hansen_projection)\n",
    "        Coconut = Coconut.updateMask(Coconut.gte(1))       # Reduce resolution based on majority interpolation (45% majority for 20 m resolution)\n",
    "        Coconut = Coconut.where(Coconut.gt(0), Coconut_Global_reclass_value)\n",
    "        Coconut = Coconut.where(ManagedForests.gt(0).And(Coconut.gt(0)), Coconut.multiply(-1))\n",
    "    \n",
    "\n",
    "    \"\"\"\n",
    "    # Load Rubber dataset\n",
    "    Rubber = Rubber.updateMask(TC_mask)\n",
    "    Rubber = Rubber.setDefaultProjection(Rubber_projection)\n",
    "    Rubber = Rubber.reduceResolution(reducer = ee.Reducer.count(), bestEffort = True, maxPixels = 3*3).reproject(Hansen_projection)\n",
    "    Rubber = Rubber.updateMask(Rubber.gt(4))       # Reduce resolution based on majority interpolation \n",
    "    Rubber = Rubber.where(Rubber.gt(0), Rubber_reclass_values)\n",
    "    Rubber = Rubber.where(ManagedForests.gt(0).And(Rubber.gt(0)), Rubber.multiply(-1))\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    ## Analyze forest fire\n",
    "    Forest_fire = Forest_fire.updateMask(TC_mask)\n",
    "    Forest_fire = Forest_fire.updateMask(Forest_fire.gt(1))\n",
    "    Forest_fire = Forest_fire.where(ManagedForests.gt(0).And(Forest_fire.gt(1)), Forest_fire.multiply(-1))\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## Analyze dominant forest loss driver\n",
    "    #dominant_driver = dominant_driver.clip(geometry)#$$\n",
    "    dominant_driver = dominant_driver.updateMask(TC_mask)\n",
    "    dominant_driver = dominant_driver.updateMask(dominant_driver.gt(0))\n",
    "    # Convert dominant_driver to negative where it does not overlap with ForestManagement\n",
    "    dominant_driver = dominant_driver.where(ManagedForests.gt(0).And(dominant_driver.gt(1)), dominant_driver.multiply(-1))\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    ####################\n",
    "    # 5. Analyse attribution (priority-based)\n",
    "    ####################\n",
    "    \n",
    "    # Convert Hansen_loss to Hansen_loss_attribution\n",
    "    Hansen_loss_attribution = Hansen_loss.rename('classification')\n",
    "    \n",
    "    # Attributing forest loss to OilPalm in Indonesia\n",
    "    if Admin_boundary == 'Indonesia':\n",
    "        # Since we are not considering Oil Palm values in Du et al plantation for Indonesia\n",
    "        FL_att_to_old_plantations_species = FL_att_to_old_plantations_species.updateMask(FL_att_to_old_plantations_species.neq(Oilpalm_Plantation_reclass_value))\n",
    "        FL_att_to_new_plantations_species = FL_att_to_new_plantations_species.where(FL_att_to_new_plantations_species.eq(Oilpalm_Plantation_reclass_value), 1)\n",
    "    \n",
    "        # Oil Palm Indonesia data is preferred over MapBiomas (directly converted in the reclassed function) or Forest plantation, thus all Oil Palm-MapBiomas pixels are converted to 1 (to include them in attribution pool)\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(OilPalm_Indonesia_Image.gt(0)), Hansen_loss_attribution.multiply(Oilpalm_Indonesia_reclass_value))\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(OilPalm_Indonesia_pre2000_Image.gt(0)), Hansen_loss_attribution.multiply(-Oilpalm_Indonesia_reclass_value))\n",
    "    \n",
    "    # Attributing forest loss to Maize-China\n",
    "    if Admin_boundary == 'China':\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(20)).And(Maize_China_four), Maize_China_four)\n",
    "\n",
    "    # Attributing forest loss to Soybean\n",
    "    if Admin_boundary in MapBiomas_dataset_version_country:\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(22)).And(Soy_four.gt(0)), Soy_four)\n",
    "\n",
    "    # Attributing forest loss to Sugarcane in Brazil\n",
    "    if Admin_boundary == 'Brazil':\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(19)).And(Sugarcane_Brazil), Sugarcane_Brazil)\n",
    "    \n",
    "    # Attributing forest loss to OilPalm-Malaysia\n",
    "    if Admin_boundary == 'Malaysia':\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(OilPalm_Malaysia_four), OilPalm_Malaysia_four)\n",
    "\n",
    "    # Attributing forest loss to Cocoa\n",
    "    if Admin_boundary in ['Côte d\\'Ivoire', 'Ghana']: \n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(21)).And(Cocoa_map_CDG), Cocoa_map_CDG)\n",
    "\n",
    "    # Attributing MapBiomas 'commodities' specifically (Match with commodities in the Supplementary data)\n",
    "    if Admin_boundary in MapBiomas_dataset_version_country:\n",
    "        if Admin_boundary == 'Bolivia':\n",
    "            MapBiomas_endyear = 21 #For 2021\n",
    "        else:\n",
    "            MapBiomas_endyear = 22 #For 2022\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(MapBiomas_endyear)).And(MapBiomas.gte(3221)).And(MapBiomas.lte(3999)), MapBiomas) \\\n",
    "                                                         .where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(MapBiomas_endyear)).And(MapBiomas.eq(4000)), MapBiomas) \\\n",
    "                                                         .where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(MapBiomas_endyear)).And(MapBiomas.gte(6001)).And(MapBiomas.lte(6999)), MapBiomas)\n",
    "    \n",
    "    # Attributing forest loss to Rice\n",
    "    if Admin_boundary in Rice_dataset_version_country:\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(19)).And(Paddy_rice), Paddy_rice)\n",
    "\n",
    "    # Attributing forest loss to Rapeseed\n",
    "    if Admin_boundary in Rapeseed_dataset_version_country:\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(19)).And(Rapeseed), Rapeseed)\n",
    "    \n",
    "    # Attributing forest loss to Rubber\n",
    "    #$$Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(22)).And(Rubber), Rubber)\n",
    "\n",
    "    if Admin_boundary not in Global_oilpalm_and_Coconut_not_in:\n",
    "        # Attributing forest loss to OilPalm (Global; except Indonesia)\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(19)).And(OilPalm_Global), OilPalm_Global)\n",
    "        \n",
    "        # Attributing forest loss to Coconut\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(20)).And(Coconut), Coconut)\n",
    "\n",
    "    # Attributing forest loss to Plantation species (old and new)\n",
    "    if Admin_boundary in Plantation_dataset_version_country:\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(FL_att_to_old_plantations_species.gt(0)), FL_att_to_old_plantations_species.multiply(-1))\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(FL_att_to_new_plantations_species.gt(0)), FL_att_to_new_plantations_species)\n",
    "        \n",
    "    # Attributing forest loss to MapBiomas landuse\n",
    "    if Admin_boundary in MapBiomas_dataset_version_country:\n",
    "        Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(MapBiomas_endyear)).And(MapBiomas), MapBiomas)\n",
    "\n",
    "    # Attributing forest loss to arable/temporary cropland\n",
    "    Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Hansen_lossyear.lte(19)).And(combinedCropland), combinedCropland)\n",
    "\n",
    "    # Attributing forest loss to Forest fire\n",
    "    Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(Forest_fire), Forest_fire)\n",
    "\n",
    "    # Attributing forest loss to Dominant drivers\n",
    "    Hansen_loss_attribution = Hansen_loss_attribution.where(Hansen_loss_attribution.eq(1).And(dominant_driver), dominant_driver)\n",
    "\n",
    "    \n",
    "    \n",
    "    ## Miscellaneous, since it is difficult to save files with special characters\n",
    "    if Admin_boundary == 'México':\n",
    "        Admin_boundary = 'Mexico'\n",
    "    if Admin_boundary == 'Saint-Barthélemy':\n",
    "        Admin_boundary = 'Saint-Barthelemy'\n",
    "    if Admin_boundary =='Côte d\\'Ivoire':\n",
    "        Admin_boundary = 'Cote dIvoire'\n",
    "    if Admin_boundary == 'Réunion':\n",
    "        Admin_boundary = 'Reunion'\n",
    "    if Admin_boundary == 'Åland':\n",
    "        Admin_boundary = 'Aland'\n",
    "    if Admin_boundary == 'São Tomé and Príncipe':\n",
    "        Admin_boundary = 'Sao Tome and Principe'\n",
    "    if Admin_boundary == 'Curaçao':\n",
    "        Admin_boundary = 'Curacao'\n",
    "    \n",
    "    \n",
    "    ####################\n",
    "    # 6. Analyse above and below ground biomass\n",
    "    ####################\n",
    "\n",
    "    if simulation in ['All', 'AGB', 'BGB']:\n",
    "        AGB = AGB.updateMask(TC_mask)\n",
    "        #$$AGB = AGB.clip(geometry)\n",
    "        AGB = AGB.multiply(44/12).multiply(0.47)  # Convert AGB to C, and then to CO2\n",
    "        \n",
    "        # Ecoregion for BGB calculation\n",
    "        ecoRegions = ee.FeatureCollection(Ecoregion_dataset_version)\n",
    "        ecoregionImage = ecoRegions.reduceToImage(\n",
    "            properties=['BIOME_NUM'], \n",
    "            reducer=ee.Reducer.sum(),\n",
    "        )\n",
    "        \n",
    "        ecoregionImage = ecoregionImage.updateMask(TC_mask)\n",
    "        #$$ecoregionImage = ecoregionImage.clip(geometry)\n",
    "        \n",
    "        ## Analyse deadwood and litter\n",
    "        elevation = elevation.updateMask(TC_mask)\n",
    "        precipitation = precipitation.updateMask(TC_mask)\n",
    "            \n",
    "        # Define conditions\n",
    "        tropical = ecoregionImage.eq(1).Or(ecoregionImage.eq(2)).Or(ecoregionImage.eq(3)).Or(ecoregionImage.eq(7))\n",
    "        lowElevation = elevation.lt(2000)\n",
    "        mediumPrecipitation = precipitation.gt(1000).And(precipitation.lte(1600))\n",
    "        highPrecipitation = precipitation.gt(1600)\n",
    "        highElevation = elevation.gt(2000)\n",
    "        temperateBoreal = tropical.Not()\n",
    "            \n",
    "        # Apply conditions\n",
    "        condition1 = tropical.And(lowElevation).And(precipitation.lt(1000))\n",
    "        DeadWood_scale1 = ee.Image.constant(0.02)\n",
    "        Litter_scale1 = ee.Image.constant(0.04)\n",
    "            \n",
    "        condition2 = tropical.And(lowElevation).And(mediumPrecipitation)\n",
    "        DeadWood_scale2 = ee.Image.constant(0.01)\n",
    "        Litter_scale2 = ee.Image.constant(0.01)\n",
    "            \n",
    "        condition3 = tropical.And(lowElevation).And(highPrecipitation)\n",
    "        DeadWood_scale3 = ee.Image.constant(0.06)\n",
    "        Litter_scale3 = ee.Image.constant(0.01)\n",
    "            \n",
    "        condition4 = tropical.And(highElevation)\n",
    "        DeadWood_scale4 = ee.Image.constant(0.07)\n",
    "        Litter_scale4 = ee.Image.constant(0.01)\n",
    "            \n",
    "        condition5 = temperateBoreal\n",
    "        DeadWood_scale5 = ee.Image.constant(0.08)\n",
    "        Litter_scale5 = ee.Image.constant(0.04)\n",
    "            \n",
    "        Deadwood = AGB.where(condition1, AGB.multiply(DeadWood_scale1)) \\\n",
    "                        .where(condition2, AGB.multiply(DeadWood_scale2)) \\\n",
    "                        .where(condition3, AGB.multiply(DeadWood_scale3)) \\\n",
    "                        .where(condition4, AGB.multiply(DeadWood_scale4)) \\\n",
    "                        .where(condition5, AGB.multiply(DeadWood_scale5))\n",
    "            \n",
    "        Litter = AGB.where(condition1, AGB.multiply(Litter_scale1)) \\\n",
    "                    .where(condition2, AGB.multiply(Litter_scale2)) \\\n",
    "                    .where(condition3, AGB.multiply(Litter_scale3)) \\\n",
    "                    .where(condition4, AGB.multiply(Litter_scale4)) \\\n",
    "                    .where(condition5, AGB.multiply(Litter_scale5))\n",
    "        \n",
    "        AGB_total = AGB.add(Deadwood).add(Litter)\n",
    "    \n",
    "        if simulation in ['All', 'BGB']:       \n",
    "            # Define the BGB/AGB multiplication factors based on ecoregion value and AGB (Scale values define thje ration converted from Mg/ha to Mg/px)\n",
    "            Root = ee.Image(\"projects/lu-chandrakant/assets/Root_Shoot_ratio/AROOT\");\n",
    "            Shoot = ee.Image(\"projects/lu-chandrakant/assets/Root_Shoot_ratio/ASHOOT\");\n",
    "                \n",
    "            BGB_AGB = Root.divide(Shoot).unmask(0);\n",
    "            \n",
    "            # Apply conditional statements using ee.Image.where()\n",
    "            BGB = AGB.multiply(BGB_AGB);\n",
    "\n",
    "\n",
    "    if simulation in ['All', 'PEATLAND', 'SOC 0-30', 'SOC 30-100']:       \n",
    "        # Clip the image to a specified geometry and apply a mask (TC_mask)\n",
    "        Peatland = Peatland.updateMask(TC_mask)\n",
    "        # Adding values from Hansen_loss_attribution to Peatland pixels\n",
    "        Peatland = Hansen_loss_attribution.updateMask(Peatland) \n",
    "        \n",
    "    \n",
    "    ## Analyse Soil Organic Carbon Stock\n",
    "    if simulation in ['All', 'SOC 0-30', 'SOC 30-100']:   \n",
    "        # Calculate SOC stock for different depth ranges\n",
    "        pixelarea = ee.Image.pixelArea().reproject(Hansen_projection).divide(10000)\n",
    "\n",
    "    if simulation in ['All', 'SOC 0-30']:   \n",
    "        # Combine SOC layers for different depth ranges\n",
    "        SOC_0_30 = SOC_0_30.updateMask(TC_mask)\n",
    "        # Converting tC/ha to tC/px\n",
    "        SOC_0_30 = SOC_0_30.where(SOC_0_30, SOC_0_30.multiply(pixelarea))\n",
    "        SOC_0_30 = SOC_0_30.multiply(44/12)  # Convert tC/px to tCO2/px\n",
    "        SOC_0_30 = SOC_0_30.where(Peatland, SOC_0_30.multiply(0))\n",
    "\n",
    "\n",
    "    if simulation in ['All', 'SOC 30-100']:   \n",
    "        SOC_30_60_t_per_ha = calculate_SOC(30, 60)\n",
    "        SOC_60_100_t_per_ha = calculate_SOC(60, 100)\n",
    "        SOC_30_100 = (SOC_30_60_t_per_ha.add(SOC_60_100_t_per_ha)).updateMask(TC_mask)    \n",
    "        SOC_30_100 = SOC_30_100.where(SOC_30_100, SOC_30_100.multiply(pixelarea))\n",
    "        SOC_30_100 = SOC_30_100.multiply(44/12)    \n",
    "        SOC_30_100 = SOC_30_100.where(Peatland, SOC_30_100.multiply(0))\n",
    "        \n",
    "    \n",
    "        \n",
    "    \n",
    "    ####################\n",
    "    # 7. Exporting attribution as csv\n",
    "    ####################\n",
    "    \n",
    "    ## Export classified forest loss attribution as csv\n",
    "    input = Hansen_loss_attribution\n",
    "\n",
    "    def Classification_to_region_features(region):\n",
    "        def to_class_feature(class_group):\n",
    "            class_group = ee.Feature(class_group)\n",
    "            year_groups = ee.List(class_group.get('groups'))\n",
    "            areas = ee.Array(\n",
    "                year_groups.map(\n",
    "                    lambda properties: ee.Dictionary(properties).toArray().toList()\n",
    "                )\n",
    "            )\n",
    "            labels = (\n",
    "                areas.slice(1, 0, 1)\n",
    "                .project([0])\n",
    "                .toList()\n",
    "                .map(\n",
    "                    lambda year: ee.Number(year).toInt().add(2000).format('loss_%d'))\n",
    "            )\n",
    "            values = areas.slice(1, 1).project([0]).toList()\n",
    "            return (\n",
    "                ee.Feature(\n",
    "                    None,\n",
    "                    region.toDictionary(['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2'])\n",
    "                )\n",
    "                .set('Class', class_group.get('Class'))\n",
    "                .set(ee.Dictionary.fromLists(labels, values))\n",
    "            )\n",
    "\n",
    "        groups = (\n",
    "            ee.FeatureCollection(\n",
    "                ee.List(\n",
    "                    ee.Image.pixelArea().reproject(Hansen_projection).divide(1e4)\n",
    "                    .addBands(Hansen_lossyear)\n",
    "                    .addBands(input)\n",
    "                    .reduceRegion(\n",
    "                        reducer=ee.Reducer.sum().group(1, 'lossYear').group(2, 'Class'),\n",
    "                        geometry=region.geometry(),\n",
    "                        scale=Hansen_lossyear.projection().nominalScale(),\n",
    "                        maxPixels=1e13,\n",
    "                        tileScale=tileScale_value,\n",
    "                    ).get('groups')\n",
    "                ).map(\n",
    "                    # Include a size, to filter out cases where no data exist\n",
    "                    lambda group: ee.Feature(\n",
    "                        None,\n",
    "                        ee.Dictionary(group)\n",
    "                    ).set(\n",
    "                        'size',\n",
    "                        ee.List(ee.Dictionary(group).get('groups')).size()\n",
    "                    )\n",
    "                )\n",
    "            )\n",
    "            .filter(ee.Filter.gt('size', 0))\n",
    "        )\n",
    "        return ee.FeatureCollection(groups.map(to_class_feature))\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## Export AGB for each forest loss class as csv\n",
    "    def AGB_Emission_to_region_features(region):\n",
    "        def to_class_feature(class_group):\n",
    "            class_group = ee.Feature(class_group)\n",
    "            year_groups = ee.List(class_group.get('groups'))\n",
    "            areas = ee.Array(\n",
    "                year_groups.map(\n",
    "                    lambda properties: ee.Dictionary(properties).toArray().toList()\n",
    "                )\n",
    "            )\n",
    "            labels = (\n",
    "                areas.slice(1, 0, 1)\n",
    "                .project([0])\n",
    "                .toList()\n",
    "                .map(\n",
    "                    lambda year: ee.Number(year).toInt().add(2000).format('loss_%d'))\n",
    "            )\n",
    "            values = areas.slice(1, 1).project([0]).toList()\n",
    "            return (\n",
    "                ee.Feature(\n",
    "                    None,\n",
    "                    region.toDictionary(['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2'])\n",
    "                )\n",
    "                .set('Class', class_group.get('Class'))\n",
    "                .set(ee.Dictionary.fromLists(labels, values))\n",
    "            )\n",
    "\n",
    "        groups = (\n",
    "            ee.FeatureCollection(\n",
    "                ee.List(\n",
    "                    AGB_total          # AGB is called here (includes Deadwood and Litter)\n",
    "                    .addBands(Hansen_lossyear)\n",
    "                    .addBands(input)\n",
    "                    .reduceRegion(\n",
    "                        reducer=ee.Reducer.sum().group(1, 'lossYear').group(2, 'Class'),\n",
    "                        geometry=region.geometry(),\n",
    "                        scale=Hansen_lossyear.projection().nominalScale(),\n",
    "                        maxPixels=1e13,\n",
    "                        tileScale=tileScale_value,\n",
    "                    ).get('groups')\n",
    "                ).map(\n",
    "                    # Include a size, to filter out cases where no data exist\n",
    "                    lambda group: ee.Feature(\n",
    "                        None,\n",
    "                        ee.Dictionary(group)\n",
    "                    ).set(\n",
    "                        'size',\n",
    "                        ee.List(ee.Dictionary(group).get('groups')).size()\n",
    "                    )\n",
    "                )\n",
    "            )\n",
    "            .filter(ee.Filter.gt('size', 0))\n",
    "        )\n",
    "        return ee.FeatureCollection(groups.map(to_class_feature))\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## Export BGB for each forest loss class as csv\n",
    "    def BGB_Emission_to_region_features(region):\n",
    "        def to_class_feature(class_group):\n",
    "            class_group = ee.Feature(class_group)\n",
    "            year_groups = ee.List(class_group.get('groups'))\n",
    "            areas = ee.Array(\n",
    "                year_groups.map(\n",
    "                    lambda properties: ee.Dictionary(properties).toArray().toList()\n",
    "                )\n",
    "            )\n",
    "            labels = (\n",
    "                areas.slice(1, 0, 1)\n",
    "                .project([0])\n",
    "                .toList()\n",
    "                .map(\n",
    "                    lambda year: ee.Number(year).toInt().add(2000).format('loss_%d'))\n",
    "            )\n",
    "            values = areas.slice(1, 1).project([0]).toList()\n",
    "            return (\n",
    "                ee.Feature(\n",
    "                    None,\n",
    "                    region.toDictionary(['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2'])\n",
    "                )\n",
    "                .set('Class', class_group.get('Class'))\n",
    "                .set(ee.Dictionary.fromLists(labels, values))\n",
    "            )\n",
    "\n",
    "        groups = (\n",
    "            ee.FeatureCollection(\n",
    "                ee.List(\n",
    "                    BGB          # BGB is called here\n",
    "                    .addBands(Hansen_lossyear)\n",
    "                    .addBands(input)\n",
    "                    .reduceRegion(\n",
    "                        reducer=ee.Reducer.sum().group(1, 'lossYear').group(2, 'Class'),\n",
    "                        geometry=region.geometry(),\n",
    "                        scale=Hansen_lossyear.projection().nominalScale(),\n",
    "                        maxPixels=1e13,\n",
    "                        tileScale=tileScale_value,\n",
    "                    ).get('groups')\n",
    "                ).map(\n",
    "                    # Include a size, to filter out cases where no data exist\n",
    "                    lambda group: ee.Feature(\n",
    "                        None,\n",
    "                        ee.Dictionary(group)\n",
    "                    ).set(\n",
    "                        'size',\n",
    "                        ee.List(ee.Dictionary(group).get('groups')).size()\n",
    "                    )\n",
    "                )\n",
    "            )\n",
    "            .filter(ee.Filter.gt('size', 0))\n",
    "        )\n",
    "        return ee.FeatureCollection(groups.map(to_class_feature))\n",
    "\n",
    "\n",
    "    ## Export SOC 0-30m for each forest loss class as csv\n",
    "    def SOC_0_30_Emission_to_region_features(region):\n",
    "        def to_class_feature(class_group):\n",
    "            class_group = ee.Feature(class_group)\n",
    "            year_groups = ee.List(class_group.get('groups'))\n",
    "            areas = ee.Array(\n",
    "                year_groups.map(\n",
    "                    lambda properties: ee.Dictionary(properties).toArray().toList()\n",
    "                )\n",
    "            )\n",
    "            labels = (\n",
    "                areas.slice(1, 0, 1)\n",
    "                .project([0])\n",
    "                .toList()\n",
    "                .map(\n",
    "                    lambda year: ee.Number(year).toInt().add(2000).format('loss_%d'))\n",
    "            )\n",
    "            values = areas.slice(1, 1).project([0]).toList()\n",
    "            return (\n",
    "                ee.Feature(\n",
    "                    None,\n",
    "                    region.toDictionary(['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2'])\n",
    "                )\n",
    "                .set('Class', class_group.get('Class'))\n",
    "                .set(ee.Dictionary.fromLists(labels, values))\n",
    "            )\n",
    "\n",
    "        groups = (\n",
    "            ee.FeatureCollection(\n",
    "                ee.List(\n",
    "                    SOC_0_30          # SOC 0-30 m\n",
    "                    .addBands(Hansen_lossyear)\n",
    "                    .addBands(input)\n",
    "                    .reduceRegion(\n",
    "                        reducer=ee.Reducer.sum().group(1, 'lossYear').group(2, 'Class'),\n",
    "                        geometry=region.geometry(),\n",
    "                        scale=Hansen_lossyear.projection().nominalScale(),\n",
    "                        maxPixels=1e13,\n",
    "                        tileScale=tileScale_value,\n",
    "                    ).get('groups')\n",
    "                ).map(\n",
    "                    # Include a size, to filter out cases where no data exist\n",
    "                    lambda group: ee.Feature(\n",
    "                        None,\n",
    "                        ee.Dictionary(group)\n",
    "                    ).set(\n",
    "                        'size',\n",
    "                        ee.List(ee.Dictionary(group).get('groups')).size()\n",
    "                    )\n",
    "                )\n",
    "            )\n",
    "            .filter(ee.Filter.gt('size', 0))\n",
    "        )\n",
    "        return ee.FeatureCollection(groups.map(to_class_feature))\n",
    "\n",
    "\n",
    "    ## Export SOC 30-100m for each forest loss class as csv\n",
    "    def SOC_30_100_Emission_to_region_features(region):\n",
    "        def to_class_feature(class_group):\n",
    "            class_group = ee.Feature(class_group)\n",
    "            year_groups = ee.List(class_group.get('groups'))\n",
    "            areas = ee.Array(\n",
    "                year_groups.map(\n",
    "                    lambda properties: ee.Dictionary(properties).toArray().toList()\n",
    "                )\n",
    "            )\n",
    "            labels = (\n",
    "                areas.slice(1, 0, 1)\n",
    "                .project([0])\n",
    "                .toList()\n",
    "                .map(\n",
    "                    lambda year: ee.Number(year).toInt().add(2000).format('loss_%d'))\n",
    "            )\n",
    "            values = areas.slice(1, 1).project([0]).toList()\n",
    "            return (\n",
    "                ee.Feature(\n",
    "                    None,\n",
    "                    region.toDictionary(['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2'])\n",
    "                )\n",
    "                .set('Class', class_group.get('Class'))\n",
    "                .set(ee.Dictionary.fromLists(labels, values))\n",
    "            )\n",
    "\n",
    "        groups = (\n",
    "            ee.FeatureCollection(\n",
    "                ee.List(\n",
    "                    SOC_30_100          # SOC 30-100 m\n",
    "                    .addBands(Hansen_lossyear)\n",
    "                    .addBands(input)\n",
    "                    .reduceRegion(\n",
    "                        reducer=ee.Reducer.sum().group(1, 'lossYear').group(2, 'Class'),\n",
    "                        geometry=region.geometry(),\n",
    "                        scale=Hansen_lossyear.projection().nominalScale(),\n",
    "                        maxPixels=1e13,\n",
    "                        tileScale=tileScale_value,\n",
    "                    ).get('groups')\n",
    "                ).map(\n",
    "                    # Include a size, to filter out cases where no data exist\n",
    "                    lambda group: ee.Feature(\n",
    "                        None,\n",
    "                        ee.Dictionary(group)\n",
    "                    ).set(\n",
    "                        'size',\n",
    "                        ee.List(ee.Dictionary(group).get('groups')).size()\n",
    "                    )\n",
    "                )\n",
    "            )\n",
    "            .filter(ee.Filter.gt('size', 0))\n",
    "        )\n",
    "        return ee.FeatureCollection(groups.map(to_class_feature))\n",
    "\n",
    "\n",
    "    def Classification_to_Peatland_regions(region):\n",
    "        def to_class_feature(class_group):\n",
    "            class_group = ee.Feature(class_group)\n",
    "            year_groups = ee.List(class_group.get('groups'))\n",
    "            areas = ee.Array(\n",
    "                year_groups.map(\n",
    "                    lambda properties: ee.Dictionary(properties).toArray().toList()\n",
    "                )\n",
    "            )\n",
    "            labels = (\n",
    "                areas.slice(1, 0, 1)\n",
    "                .project([0])\n",
    "                .toList()\n",
    "                .map(\n",
    "                    lambda year: ee.Number(year).toInt().add(2000).format('loss_%d'))\n",
    "            )\n",
    "            values = areas.slice(1, 1).project([0]).toList()\n",
    "            return (\n",
    "                ee.Feature(\n",
    "                    None,\n",
    "                    region.toDictionary(['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2'])\n",
    "                )\n",
    "                .set('Class', class_group.get('Class'))\n",
    "                .set(ee.Dictionary.fromLists(labels, values))\n",
    "            )\n",
    "\n",
    "        groups = (\n",
    "            ee.FeatureCollection(\n",
    "                ee.List(\n",
    "                    ee.Image.pixelArea().reproject(Hansen_projection).divide(1e4)\n",
    "                    .addBands(Hansen_lossyear)\n",
    "                    .addBands(Peatland)                 #Changed data to Peatland for attribution\n",
    "                    .reduceRegion(\n",
    "                        reducer=ee.Reducer.sum().group(1, 'lossYear').group(2, 'Class'),\n",
    "                        geometry=region.geometry(),\n",
    "                        scale=Hansen_lossyear.projection().nominalScale(),\n",
    "                        maxPixels=1e13,\n",
    "                        tileScale=tileScale_value,\n",
    "                    ).get('groups')\n",
    "                ).map(\n",
    "                    # Include a size, to filter out cases where no data exist\n",
    "                    lambda group: ee.Feature(\n",
    "                        None,\n",
    "                        ee.Dictionary(group)\n",
    "                    ).set(\n",
    "                        'size',\n",
    "                        ee.List(ee.Dictionary(group).get('groups')).size()\n",
    "                    )\n",
    "                )\n",
    "            )\n",
    "            .filter(ee.Filter.gt('size', 0))\n",
    "        )\n",
    "        return ee.FeatureCollection(groups.map(to_class_feature))\n",
    "    \n",
    "    \n",
    "    # Only exports classification\n",
    "    def Forest_loss_attribution(region, LOSS, runID = ''):\n",
    "        if LOSS == 'CLASSIFICATION':\n",
    "            areas = region.map(Classification_to_region_features).flatten()\n",
    "            savefilename = 'CLASSIFICATION'\n",
    "        elif LOSS == 'PEATLAND':\n",
    "            areas = region.map(Classification_to_Peatland_regions).flatten()\n",
    "            savefilename = 'PEATLAND'\n",
    "        propertiesToExport = ['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2', 'Class']\n",
    "        for i in range(2001, end_year+1):\n",
    "            propertiesToExport.append('loss_' + str(i))\n",
    "\n",
    "        task = ee.batch.Export.table.toDrive(\n",
    "            collection=areas,\n",
    "            description= 'Forest_loss_to_'+savefilename+'_'+ str(Admin_boundary)+'_'+ runID + '_' +'PythonAPI_'+ str(datetime.datetime.now().isoformat()[0:19]),\n",
    "            folder='DeDuCE_v'+str(GEE_version)+sub_folder_code,\n",
    "            fileFormat='CSV',\n",
    "            selectors=propertiesToExport\n",
    "        )\n",
    "        task.start()\n",
    "        \n",
    "    # Exports AGB and BGB\n",
    "    def Forest_loss_emission_attribution(region, emission, runID = ''):\n",
    "        if emission == 'BGB':\n",
    "            areas = region.map(BGB_Emission_to_region_features).flatten()\n",
    "            savefilename = 'BGB_EMISSION'\n",
    "        elif emission == 'SOC 0-30':\n",
    "            areas = region.map(SOC_0_30_Emission_to_region_features).flatten()\n",
    "            savefilename = 'SOC_0_30'\n",
    "        elif emission == 'SOC 30-100':\n",
    "            areas = region.map(SOC_30_100_Emission_to_region_features).flatten()\n",
    "            savefilename = 'SOC_30_100'\n",
    "        elif emission == 'AGB':\n",
    "            areas = region.map(AGB_Emission_to_region_features).flatten()\n",
    "            savefilename = 'AGB_EMISSION'\n",
    "        propertiesToExport = ['CONTINENT', 'COUNTRY', 'GID_0', 'GID_1', 'GID_2', 'Class']\n",
    "        for i in range(2001, end_year+1):\n",
    "            propertiesToExport.append('loss_' + str(i))\n",
    "\n",
    "        task = ee.batch.Export.table.toDrive(\n",
    "            collection=areas,\n",
    "            description= 'Forest_loss_to_'+savefilename+'_'+ str(Admin_boundary)+'_'+ runID + '_' +'PythonAPI_'+ str(datetime.datetime.now().isoformat()[0:19]),\n",
    "            folder='DeDuCE_v'+str(GEE_version)+sub_folder_code,\n",
    "            fileFormat='CSV',\n",
    "            selectors=propertiesToExport\n",
    "        )\n",
    "        task.start()\n",
    "    \n",
    "    \n",
    "    ####################\n",
    "    # Funtions that run the export\n",
    "    ####################\n",
    "    tileScale_value = 2\n",
    "\n",
    "    if simulation == 'CLASSIFICATION':\n",
    "        Forest_loss_attribution(geometry, 'CLASSIFICATION')\n",
    "    elif simulation == 'PEATLAND':\n",
    "        Forest_loss_attribution(geometry, 'PEATLAND')\n",
    "    elif simulation == 'AGB':\n",
    "        Forest_loss_emission_attribution(geometry, 'AGB')\n",
    "    elif simulation == 'BGB':\n",
    "        Forest_loss_emission_attribution(geometry, 'BGB')\n",
    "    elif simulation == 'SOC 0-30':\n",
    "        Forest_loss_emission_attribution(geometry, 'SOC 0-30')\n",
    "    elif simulation == 'SOC 30-100':\n",
    "        Forest_loss_emission_attribution(geometry, 'SOC 30-100')\n",
    "    else:\n",
    "        Forest_loss_attribution(geometry, 'CLASSIFICATION')\n",
    "        Forest_loss_attribution(geometry, 'PEATLAND')\n",
    "        Forest_loss_emission_attribution(geometry, 'AGB')\n",
    "        Forest_loss_emission_attribution(geometry, 'BGB')\n",
    "        Forest_loss_emission_attribution(geometry, 'SOC 0-30')\n",
    "        Forest_loss_emission_attribution(geometry, 'SOC 30-100')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3cbb0cba-3c42-45a0-861b-0199ac78c125",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Boundary Uruguay : 141\n"
     ]
    }
   ],
   "source": [
    "GEE_loop('Uruguay')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "2dd44a13-1c6e-42ca-8b2b-673775b60b87",
   "metadata": {},
   "source": [
    "for Country in Countries.COUNTRY.values:\n",
    "    if Country in ['Greenland', 'Kiribati']:\n",
    "        continue\n",
    "    GEE_loop(Country)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "62748505-61dc-4e29-aa88-87a1730407b9",
   "metadata": {},
   "source": [
    "import time\n",
    "import re\n",
    "\n",
    "def check_tasks_ready():\n",
    "    global tasks, failed_tasks\n",
    "    while True:\n",
    "        tasks = ee.batch.Task.list()\n",
    "        tasks = tasks[::-1]\n",
    "        all_not_ready = all(task.status()['state'] not in ['RUNNING', 'READY'] for task in ee.batch.Task.list()[:30]) #$$ Remove [:5]\n",
    "\n",
    "        if all_not_ready:\n",
    "            print('1st iteration of all tasks are completed')\n",
    "            failed_tasks = [task for task in tasks if task.status()['state'] == 'FAILED']\n",
    "\n",
    "            for task in failed_tasks:\n",
    "                task_name = task.status()['description']\n",
    "                match = re.search(r\"Forest_loss_to_((?:[^_]+_)*[^_]+)_([^_]+(?: [^_]+)*)__\", task_name)\n",
    "                \n",
    "                if match:\n",
    "                    simulation, country = match.groups()\n",
    "                    simulation = simulation.replace('SOC_30_100', 'SOC 30-100').replace('SOC_0_30', 'SOC 0-30').replace('BGB_EMISSION', 'BGB').replace('AGB_EMISSION', 'AGB')\n",
    "                    print(f\"Country: {country}, Simulation: {simulation}\")\n",
    "                    if country in ['Russia', 'United States', 'Romania', 'Mexico', \n",
    "                                   'Thailand','Algeria','Colombia','Vietnam', 'Philippines',\n",
    "                                   'Canada','Brazil',\n",
    "                                   'Svalbard and Jan Mayen']:continue\n",
    "                    while True:\n",
    "                        #if ee.batch.Task.list()[0].status()['state'] != 'READY': #$$ Change to 'READY'\n",
    "                        if all(task.status()['state'] not in ['RUNNING', 'READY'] for task in ee.batch.Task.list()[:30]):\n",
    "                            country_corrections = {\n",
    "                                'Mexico': 'México',\n",
    "                                'Saint-Barthelemy': 'Saint-Barthélemy',\n",
    "                                'Cote dIvoire': 'Côte d\\'Ivoire',\n",
    "                                'Reunion': 'Réunion',\n",
    "                                'Aland': 'Åland',\n",
    "                                'Sao Tome and Principe': 'São Tomé and Príncipe',\n",
    "                                'Curacao': 'Curaçao'\n",
    "                            }\n",
    "                            \n",
    "                            # Update the country name using the dictionary\n",
    "                            country = country_corrections.get(country, country)                               \n",
    "                            print(country, simulation)\n",
    "                            GEE_loop(country, simulation)\n",
    "                            break  # Ensures the loop does not become infinite\n",
    "                        else:\n",
    "                            time.sleep(60)  # Waits 2 minutes before re-checking, as suggested\n",
    "                else:\n",
    "                    print('\\033[91mData not in proper format\\033[91m')\n",
    "\n",
    "            break  # Exit the outer loop since tasks are not 'READY'\n",
    "        else:\n",
    "            print(\"Not all tasks are completed. Checking again in 20 minutes...\")\n",
    "            time.sleep(1200)  # Delay for 1200 seconds or 20 minutes\n",
    "\n",
    "# Run the function\n",
    "check_tasks_ready()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e81fe2c2-1b4c-4cfb-b51f-257be611e668",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Japan CLASSIFICATION\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "for Country in ['United States']:\n",
    "    for simulation in ['CLASSIFICATION','PEATLAND','AGB','BGB','SOC 0-30','SOC 30-100']:\n",
    "        print(Country, simulation)\n",
    "        # Load the feature collection needed to simulate aggregation\n",
    "        while True:\n",
    "            if all(task.status()['state'] not in ['RUNNING', 'READY'] for task in ee.batch.Task.list()[:5]):\n",
    "                GEE_loop(Country, simulation)\n",
    "                break\n",
    "            else:\n",
    "                #print('Waiting')\n",
    "                time.sleep(20)  # Waits 2 minutes before re-checking, as suggested"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc17fb58-97a7-4cd0-935a-24ce8a912c4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "['Romania', #Tile scale = 4; GID_0; SOC's\n",
    " 'Japan','México','Thailand','Algeria','Colombia', 'Vietnam','Philippines', #Tile scale = 2; GID_0\n",
    " 'United States','Russia','Canada','Brazil'] #Tile scale = 2; GID_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad17cd06-d661-4d4e-ad75-4f5a590a1b10",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  },
  "toc-autonumbering": true,
  "toc-showcode": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}