decision_theater_GEE_queries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Queries for data processing for CI/ASU Decision Theater presentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import GEE api and json library. Note that `ee.initialize` will fail if the key for GEE has not yet been setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import json\n",
    "import io\n",
    "\n",
    "service_account = 'gef-ldmp-server@gef-ld-toolbox.iam.gserviceaccount.com'\n",
    "credentials = ee.ServiceAccountCredentials(service_account, 'dt_key.json')\n",
    "ee.Initialize(credentials)\n",
    "#ee.Initialize()\n",
    "\n",
    "# Define a geojson as a string to use for testing\n",
    "geojson_text = '''{\"type\": \"Polygon\", \"coordinates\": [ [ [ 36.898123159960932, -0.220252698724199 ], [ 37.40818204121706, 0.915787536800825 ], [ 38.613775760549743, 0.011592247301316 ], [ 38.219639352306366, -0.892603042198193 ], [ 36.898123159960932, -0.220252698724199 ] ] ] }'''\n",
    "aoi = ee.Geometry(json.loads(geojson_text))\n",
    "\n",
    "out = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get statistics on areas degraded within a polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to pull areas that are saved as properties within a feature class,\n",
    "# convert them to percentages of the total area, and return as a dictionary. Sums\n",
    "# all features together. Scaling converts to percentages if set to 100 and\n",
    "# scaling is True\n",
    "def fc_areas_to_pct_dict(fc, normalize=True, scaling=100):\n",
    "    # Note that there may be multiple features\n",
    "    ret = {}\n",
    "    for p in [feature['properties'] for feature in fc.getInfo()['features']]:\n",
    "        areas = {}\n",
    "        # If there is more than one feature, need to update ret with these values\n",
    "        for key, value in p.iteritems():\n",
    "            if key in ret:\n",
    "                ret[key] += value\n",
    "            else:\n",
    "                ret[key] = value\n",
    "    if normalize:\n",
    "        denominator = sum(ret.values())\n",
    "        ret = {key: (value/denominator)*scaling for key, value in ret.iteritems()}\n",
    "    return ret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# polygon area in hectares\n",
    "out['area'] = aoi.area().divide(10000).getInfo()\n",
    "\n",
    "# s2_02: Number of people living inside the polygon in 2015\n",
    "pop_cnt = ee.Image(\"CIESIN/GPWv4/unwpp-adjusted-population-count/2015\");\n",
    "population = pop_cnt.reduceRegion(reducer=ee.Reducer.sum(), geometry=aoi, scale=1000, maxPixels=1e9)\n",
    "out['population'] = population.getInfo()['population-count']\n",
    "\n",
    "# s2_03: Main livelihoods\n",
    "liv = ee.FeatureCollection(\"users/geflanddegradation/toolbox_datasets/livelihoodzones\")\n",
    "\n",
    "livImage = liv.filter(ee.Filter.neq('lztype_num', None)).reduceToImage(properties=['lztype_num'], reducer=ee.Reducer.first())\n",
    "\n",
    "fields = [\"Agro-Forestry\", \"Agro-Pastoral\", \"Arid\", \"Crops - Floodzone\", \"Crops - Irrigated\", \"Crops - Rainfed\", \"Fishery\", \"Forest-Based\", \"National Park\", \"Other\", \"Pastoral\", \"Urban\"]\n",
    "# multiply pixel area by the area which experienced each of the five transitions --> output: area in ha\n",
    "livelihoodareas = livImage.eq([1,2,3,4,5,6,7,8,9,10,11,12]).rename(fields).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum(),30)\n",
    "out['livelihoods'] = fc_areas_to_pct_dict(livelihoodareas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# s3_01: SDG 15.3.1 degradation classes \n",
    "\n",
    "te_sdgi = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_sdg1531_gpg_globe_2001_2015_modis\")\n",
    "sdg_areas = te_sdgi.eq([-32768,-1,0,1]).rename([\"nodata\", \"degraded\", \"stable\", \"improving\"]).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['area_sdg'] = fc_areas_to_pct_dict(sdg_areas)\n",
    "\n",
    "# s3_02: Productivity degradation classes\n",
    "te_prod = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_lp7cl_globe_2001_2015_modis\").remap([-32768,1,2,3,4,5,6,7],[-32768,-1,-1,0,0,0,1,1])\n",
    "prod_areas = te_prod.eq([-32768,-1,0,1]).rename([\"nodata\", \"decline\", \"stable\", \"improvement\"]).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['area_prod'] = fc_areas_to_pct_dict(prod_areas)\n",
    "\n",
    "# s3_03: Land cover degradation classes\n",
    "te_land = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_lc_traj_globe_2001-2001_to_2015\")\n",
    "lc_areas = te_land.select(\"lc_dg\").eq([-1,0,1]).rename([\"degradation\", \"stable\", \"improvement\"]).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['area_lc'] = fc_areas_to_pct_dict(lc_areas)\n",
    "\n",
    "# s3_04: soc degradation classes\n",
    "te_socc_deg = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_soc_globe_2001-2015_deg\").select(\"soc_deg\")\n",
    "soc_areas = te_socc_deg.eq([-32768,-1,0,1]).rename([\"no data\", \"degradation\", \"stable\", \"improvement\"]).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['area_soc'] = fc_areas_to_pct_dict(soc_areas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# s3_05: compute land cover classes for 2001 and 2015, and the transitions which occured\n",
    "te_land = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_lc_traj_globe_2001-2001_to_2015\").select(\"lc_tr\")\n",
    "te_prod = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_lp7cl_globe_2001_2015_modis\").remap([-32768,1,2,3,4,5,6,7],[-32768,-1,-1,0,0,0,1,1])\n",
    "\n",
    "fields = [\"nodata\", \"decline\", \"stable\", \"improvement\"]\n",
    "\n",
    "prod_forests = te_prod.updateMask(te_land.eq(11)).eq([-32768,-1,0,1]).rename(fields).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['prod_forests'] = fc_areas_to_pct_dict(prod_forests)\n",
    "\n",
    "prod_grasslands = te_prod.updateMask(te_land.eq(22)).eq([-32768,-1,0,1]).rename(fields).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['prod_grasslands'] = fc_areas_to_pct_dict(prod_grasslands)\n",
    "\n",
    "prod_agriculture = te_prod.updateMask(te_land.eq(33)).eq([-32768,-1,0,1]).rename(fields).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['prod_agriculture'] = fc_areas_to_pct_dict(prod_forests)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# s3_06: compute land cover classes for 2001 and 2015, and the transitions which occured\n",
    "te_land = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_lc_traj_globe_2001-2001_to_2015\")\n",
    "\n",
    "# field names for annual land covers\n",
    "fields = [\"forest\", \"grassland\", \"agriculture\", \"wetlands\", \"artificial\", \"other land-bare\", \"water\"]\n",
    "\n",
    "# multiply pixel area by the area which experienced each of the lc classes --> output: area in ha\n",
    "lc_baseline = te_land.select(\"lc_bl\").eq([1,2,3,4,5,6,7]).rename(fields).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "# multiply pixel area by the area which experienced each of the lc classes --> output: area in ha\n",
    "lc_target = te_land.select(\"lc_tg\").eq([1,2,3,4,5,6,7]).rename(fields).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "\n",
    "out['lc_2001'] = fc_areas_to_pct_dict(lc_baseline)\n",
    "out['lc_2015'] = fc_areas_to_pct_dict(lc_target)\n",
    "\n",
    "# field names for land cover transitions between 2001-2015\n",
    "fields = [\"for-for\", \"for-gra\", \"for-agr\", \"for-wet\", \"for-art\", \"for-oth\", \"for-wat\",\n",
    "          \"gra-for\", \"gra-gra\", \"gra-agr\", \"gra-wet\", \"gra-art\", \"gra-oth\", \"gra-wat\",\n",
    "          \"agr-for\", \"agr-gra\", \"agr-agr\", \"agr-wet\", \"agr-art\", \"agr-oth\", \"agr-wat\",\n",
    "          \"wet-for\", \"wet-gra\", \"wet-agr\", \"wet-wet\", \"wet-art\", \"wet-oth\", \"wet-wat\",\n",
    "          \"art-for\", \"art-gra\", \"art-agr\", \"art-wet\", \"art-art\", \"art-oth\", \"art-wat\",\n",
    "          \"oth-for\", \"oth-gra\", \"oth-agr\", \"oth-wet\", \"oth-art\", \"oth-oth\", \"oth-wat\",\n",
    "          \"wat-for\", \"wat-gra\", \"wat-agr\", \"wat-wet\", \"wat-art\", \"wat-oth\", \"wat-wat\"]\n",
    "\n",
    "# multiply pixel area by the area which experienced each of the lc transition classes --> output: area in ha\n",
    "lc_transitions = te_land.select(\"lc_tr\").eq([11,12,13,14,15,16,17,21,22,23,24,25,26,26,31,32,33,34,35,36,37,41,42,43,44,45,46,47,51,52,53,54,55,56,57,\n",
    "              61,62,63,64,65,66,67,71,72,73,74,75,76,77]).rename(fields).multiply(ee.Image.pixelArea().divide(10000)).reduceRegions(aoi, ee.Reducer.sum())\n",
    "out['lc_transition'] = fc_areas_to_pct_dict(lc_transitions, normalize=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# s3_07: percent change in soc stocks between 2001-2015\n",
    "soc_pch_img = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_soc_globe_2001-2015_deg\").select(\"soc_pch\")\n",
    "\n",
    "# compute statistics for region\n",
    "soc_pch = soc_pch_img.reduceRegion(reducer=ee.Reducer.mean(), geometry=aoi, scale=250, maxPixels=1e9)\n",
    "# Multiple by 100 to convert to a percentage\n",
    "out['soc_change_percent'] = soc_pch.getInfo()['soc_pch'] * 100\n",
    "  \n",
    "# s3_08: change in soc stocks in tons of co2 eq between 2001-2015\n",
    "soc_an_img = ee.Image(\"users/geflanddegradation/global_ld_analysis/r20180821_soc_globe_2001-2015_annual_soc\")\n",
    "\n",
    "# compute change in SOC between 2001 and 2015 converted to co2 eq\n",
    "soc_chg_an = (soc_an_img.select('y2015').subtract(soc_an_img.select('y2001'))).multiply(250*250/10000).multiply(3.67)\n",
    "# compute statistics for the region\n",
    "soc_chg_tons_co2e = soc_chg_an.reduceRegion(reducer=ee.Reducer.sum(), geometry=aoi, scale=250, maxPixels=1e9)\n",
    "out['soc_change_tons_co2e'] = soc_chg_tons_co2e.getInfo()['y2015']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndentationError",
     "evalue": "expected an indented block (<ipython-input-8-7847e7045098>, line 3)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-8-7847e7045098>\"\u001b[1;36m, line \u001b[1;32m3\u001b[0m\n\u001b[1;33m    print out\u001b[0m\n\u001b[1;37m        ^\u001b[0m\n\u001b[1;31mIndentationError\u001b[0m\u001b[1;31m:\u001b[0m expected an indented block\n"
     ]
    }
   ],
   "source": [
    "with io.open('out.json', 'w', encoding='utf-8') as f:\n",
    "\n",
    "print out\n",
    "  f.write(json.dumps(d, ensure_ascii=False, indent=4))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}