examples/notebooks/02_demo_using_GEE_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a174eb6f",
   "metadata": {},
   "source": [
    "<a href=\"https://githubtocolab.com/kaust-halo/geeet/blob/master/examples/notebooks/02_demo_using_GEE_data.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6025a8b6",
   "metadata": {},
   "source": [
    "# Demo using GEE data \n",
    "\n",
    "This notebook is a more realistic example (but still a toy example) demonstrating the use of the PT-JPL model with real datasets from GEE. We will use NDVI data from the [MODIS/MCD43A4_006_NDVI](https://developers.google.com/earth-engine/datasets/catalog/MODIS_MCD43A4_006_NDVI) dataset and meteorological data (including net radiation) from the [ECMWF/ERA5_LAND/HOURLY](https://developers.google.com/earth-engine/datasets/catalog/ECMWF_ERA5_LAND_HOURLY#bands) dataset. \n",
    "\n",
    "The notebook contains some cells intended to be a template for exporting the data as an asset. Other cells contain interactive visualizations of the data using `geemap`. To run these cells, you will need to install [geemap](https://geemap.org/installation/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5045019e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the following line to install geeet if needed:\n",
    "#!pip install git+https://github.com/kaust-halo/geeet\n",
    "# Use the following line to install geemap if needed:\n",
    "#!pip install geemap\n",
    "import ee\n",
    "#ee.Authenticate() # Uncomment if using Google Colab or first time using EE on this device. \n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09d5a37d",
   "metadata": {},
   "source": [
    "## Data setup \n",
    "\n",
    "This section only defines the datasets, no actual GEE computations are done yet. Alternatively, skip the next two cells and load the image as an asset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f103567",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Location and time\n",
    "date_start = ee.Date('2019-06-15')\n",
    "date_end = date_start.advance(1, 'day')\n",
    "lonlat=[38.25, 30.25] \n",
    "pt1 = ee.Geometry.Point(lonlat[0]-0.25, lonlat[1]-0.25)\n",
    "pt2 = ee.Geometry.Point(lonlat[0]+0.25, lonlat[1]-0.25)\n",
    "pt3 = ee.Geometry.Point(lonlat[0]+0.25, lonlat[1]+0.25)\n",
    "pt4 = ee.Geometry.Point(lonlat[0]-0.25, lonlat[1]+0.25)\n",
    "region = ee.Geometry.Polygon([pt1,pt2,pt3,pt4])\n",
    "\n",
    "## Meteorological data: ECMWF ERA5 hourly reanalysis at 0.25° resolution. \n",
    "# For PT-JPL, we need: air temperature, relative humidity, surface pressure, and net radiation\n",
    "# https://developers.google.com/earth-engine/datasets/catalog/ECMWF_ERA5_LAND_HOURLY#bands)\n",
    "# We can define relative humidity from dewpoint temperature, air temperature, and surface pressure\n",
    "# see: https://www.ecmwf.int/sites/default/files/elibrary/2016/17117-part-iv-physical-processes.pdf\n",
    "# equation 7.91\n",
    "# net radiation = net solar radiation + net thermal radiation\n",
    "bands_to_keep = ee.List(['temperature_2m', 'dewpoint_temperature_2m', 'surface_pressure', \n",
    "    'surface_net_solar_radiation_hourly', 'surface_net_thermal_radiation_hourly'])\n",
    "bands_to_rename = ee.List(['air_temperature', 'dewpoint_temperature', 'surface_pressure', \n",
    "    'short_net_radiation_hourly', 'long_net_radiation_hourly']) \n",
    "\n",
    "Meteo_collection = (ee.ImageCollection('ECMWF/ERA5_LAND/HOURLY') \n",
    "    .filterDate(date_start, date_end) \n",
    "    .filterBounds(region)\n",
    "    .select(bands_to_keep, bands_to_rename)  \n",
    ")\n",
    "def prepare_ECMWF(img):\n",
    "    net = img.select('short_net_radiation_hourly').add(img.select('long_net_radiation_hourly'))\n",
    "    Rn = net.divide(3600).rename('net_radiation') # from J/m2 to W/m2 \n",
    "    T = img.select('air_temperature')  # in Kelvin\n",
    "    DT = img.select('dewpoint_temperature') # in Kelvin\n",
    "    P = img.select('surface_pressure') # in Pascals\n",
    "    from geeet.meteo import relative_humidity # eq 7.91 in ECMWF (2016)\n",
    "    RH = relative_humidity(T, DT)  \n",
    "    img = img.addBands(Rn).addBands(RH)\n",
    "    return(img)\n",
    "\n",
    "Meteo_collection = Meteo_collection.map(prepare_ECMWF)\n",
    "# We will work with the 10AM data as an example:\n",
    "Meteo_img = Meteo_collection.filter(ee.Filter.eq('hour',10)).first()\n",
    "\n",
    "## MODIS NDVI\n",
    "# We first start with an annual time series of NDVI\n",
    "# so we can calculate FAPAR and FAPAR_max\n",
    "NDVI_one_year = (ee.ImageCollection('MODIS/MCD43A4_006_NDVI') \n",
    "                  .filterDate(date_start.advance(-1,'year'), date_start.advance(1,'year')) \n",
    "                  .filterBounds(region)\n",
    "                  .select(\"NDVI\")\n",
    ")\n",
    "# this function calculates fapar from NDVI and adds it as a new band\n",
    "from geeet.vegetation import add_fapar\n",
    "# we can now calculate fapar_max for this period:\n",
    "fapar_max = NDVI_one_year.map(add_fapar).select('fapar').reduce(ee.Reducer.max())\n",
    "# Actual image we will use for NDVI:\n",
    "NDVI = NDVI_one_year.filterDate(date_start, date_end).first().select('NDVI')\n",
    "\n",
    "## Gather all input data into a single image:\n",
    "inputs_image = Meteo_img.addBands(NDVI).addBands(fapar_max)\n",
    "# PT-JPL requires day of year and time information:\n",
    "d = ee.Date(inputs_image.date())\n",
    "doy = d.getRelative('day', 'year').add(1)   # offset of 1 because getRelative gives 0 for doy=1\n",
    "time = d.getRelative('hour', 'day').add(3)  # offset of 3 hours to obtain local time in this region.\n",
    "inputs_image = inputs_image.set({'doy':doy,'time':time})\n",
    "inputs_image = inputs_image.clip(region)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f968b77",
   "metadata": {},
   "source": [
    "Use the next cell as a template to export the image as an asset (or skip this cell). This will perform the computations defined above to generate the input data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f838a31b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Export the input data as an earth engine asset:\n",
    "task = ee.batch.Export.image.toAsset(image=inputs_image,\n",
    "                                     description='Inputs_export_demo',\n",
    "                                     assetId='projects/geeet-public/assets/ptjpl_sample_inputs_coarse',  # Change to your own project\n",
    "                                     region=region,\n",
    "                                     scale=1000)\n",
    "#task.start()  # uncomment to submit the task."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caf6a978",
   "metadata": {},
   "source": [
    "If you skipped the last two cells, use the following to load the inputs asset:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08ca2cf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs_image = ee.Image('projects/geeet-public/assets/ptjpl_sample_inputs_coarse')\n",
    "lonlat=[38.25, 30.25] \n",
    "pt1 = ee.Geometry.Point(lonlat[0]-0.25, lonlat[1]-0.25)\n",
    "pt2 = ee.Geometry.Point(lonlat[0]+0.25, lonlat[1]-0.25)\n",
    "pt3 = ee.Geometry.Point(lonlat[0]+0.25, lonlat[1]+0.25)\n",
    "pt4 = ee.Geometry.Point(lonlat[0]-0.25, lonlat[1]+0.25)\n",
    "region = ee.Geometry.Polygon([pt1,pt2,pt3,pt4])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49c2cd68",
   "metadata": {},
   "source": [
    "Use the following cell to display an interactive map with the data inputs. Alternatively, you can [click here](https://code.earthengine.google.com/?scriptPath=users%2Flopezvoliver%2Fgeeet%3Aptjpl_sample_inputs_coarse) to open the code editor sample script. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2034bafe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optionally display the input data as an interactive map\n",
    "# (requires geeemap)\n",
    "# Alternatively, use this sample code editor script to visualize it\n",
    "# https://code.earthengine.google.com/?scriptPath=users%2Flopezvoliver%2Fgeeet%3Aptjpl_sample_inputs_coarse\n",
    "Map = None\n",
    "try:\n",
    "    import geemap\n",
    "    import geemap.colormaps as cm\n",
    "    ndvi_pal = cm.palettes.ndvi\n",
    "    pal = cm.palettes.YlOrRd\n",
    "    Map = geemap.Map(center=[lonlat[1], lonlat[0]], zoom=12)\n",
    "    Map.addLayer(inputs_image.select('relative_humidity'), {'min':0, 'max':100, 'palette':[\"440154\",\"3a528b\",\"20908d\",\"5dc962\",\"fde725\"]}, 'Relative humidity')\n",
    "    Map.addLayer(inputs_image.select('air_temperature').subtract(273.15), {'min':20, 'max':50, 'palette':pal}, 'Air temperature (C)')\n",
    "    Map.addLayer(inputs_image.select('surface_pressure').divide(1000), {'min':90, 'max':110, 'palette':pal}, 'Surface pressure (KPa)')\n",
    "    Map.addLayer(inputs_image.select('net_radiation'), {'min':0, 'max':500, 'palette':pal}, 'Net radiation (W/m2)')\n",
    "    Map.addLayer(inputs_image.select('fapar_max'), {'min':0, 'max':1, 'palette':pal}, 'FAPAR max')\n",
    "    Map.addLayer(inputs_image.select('NDVI'), {'min':0, 'max':1, 'palette':ndvi_pal}, 'NDVI')\n",
    "except:\n",
    "    pass\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08564c13",
   "metadata": {},
   "source": [
    "## Define the PT-JPL model\n",
    "\n",
    "Once the inputs image is prepared, the model can be defined simply with one line of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8768b5bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geeet.ptjpl import ptjpl_arid\n",
    "ET = ptjpl_arid(inputs_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "188606bc",
   "metadata": {},
   "source": [
    "`ET` is now an `ee.Image` object containing all of the input bands as well as the output bands. However, note that none of the model computations have run yet - not even the computations used to define the input images! (Unless of course you exported the inputs image as an asset). This is all done upon requesting data. This can be either by exporting the image as an asset (e.g. see template above), to drive, or by retrieving data using the `.getInfo()` method (e.g. see `01_geeet.ipynb`). The next cell adds the output bands (e.g., sensible and latent heat flux) to the map created above (if it exists - requires `geemap`). Alternatively, export the output image as an asset, or simply use this script to visualize the already processed output: https://code.earthengine.google.com/?scriptPath=users%2Flopezvoliver%2Fgeeet%3Aptjpl_sample_outputs_coarse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dddc473e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adds the output data to the map above. Requires \"Map\" to be defined above as a geemap.Map object. \n",
    "try:\n",
    "    import geemap.colormaps as cm\n",
    "    pal_viridis = cm.palettes.viridis\n",
    "    Map.addLayer(ET.select(['H']), {'min':0, 'max':500, 'palette':pal_viridis}, 'Sensible heat flux (W/m2)', False)\n",
    "    Map.addLayer(ET.select(['G']), {'min':0, 'max':500, 'palette':pal_viridis}, 'Ground heat flux (W/m2)', False)\n",
    "    Map.addLayer(ET.select(['Rn']), {'min':0, 'max':500, 'palette':pal_viridis}, 'Net radiation', False)\n",
    "    Map.addLayer(ET.select(['LEs']), {'min':0, 'max':500, 'palette':pal_viridis}, 'LE soil (W/m2)', False)\n",
    "    Map.addLayer(ET.select(['LEc']), {'min':0, 'max':500, 'palette':pal_viridis}, 'LE canopy (W/m2)', False)\n",
    "    Map.addLayer(ET.select(['LEi']), {'min':0, 'max':500, 'palette':pal_viridis}, 'LE interception (W/m2)', False)\n",
    "    Map.addLayer(ET.select(['LE']), {'min':0, 'max':500, 'palette':pal_viridis}, 'Total LE (W/m2)')\n",
    "except:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb49c32a",
   "metadata": {},
   "source": [
    "Use this next cell to export the output data as an asset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b033e595",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Export the output data as an earth engine asset:\n",
    "task = ee.batch.Export.image.toAsset(image=ET,\n",
    "                                     description='Outputs_export_demo',\n",
    "                                     assetId='projects/geeet-public/assets/ptjpl_sample_outputs_coarse',  # Change to your own project\n",
    "                                     region=region,\n",
    "                                     scale=1000)\n",
    "#task.start()  # uncomment to submit the task."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2f3e038",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "You can print the references for the PT-JPL model using `geeet.ptjpl.cite_all()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b2f3513",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geeet.ptjpl import cite_all\n",
    "cite_all()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "563104533abed62ee06bb251bf36028864c8de68f4f10a11edc7fddcf9524038"
  },
  "kernelspec": {
   "display_name": "Python 3.9.6 64-bit ('geeet': conda)",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}