examples/notebooks/05_xarray_landsat_era5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://githubtocolab.com/kaust-halo/geeet/blob/master/examples/notebooks/05_xarray_landsat_era5.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\"/></a>\n",
    "\n",
    "# Xarray Landsat + ERA5 ET model run (TSEB)\n",
    "\n",
    "---\n",
    "\n",
    "Here we show how to run an ET model using input data from a GeoTiff image. The image was prepared as a Cloud-Optimized Geotiff (COG), so you can read it directly from the [url](https://raw.githubusercontent.com/kaust-halo/geeet/main/examples/data/landsat_era5_sample.tif), or you can download it. \n",
    "\n",
    "The data contains Landsat-derived data (`NDVI`, `albedo`, `radiometric_temperature`), as well as climate reanalysis (ERA5) data (`surface_pressure`, `air_temperature`, `wind_speed`, `solar_radiation`, and `thermal_radiation`). The data was exported for the same image (`LC09_172039_20230421`) and region shown at the end of notebook `03_eepredefined_landsat_era5.ipynb` (or you can use [this code editor script](https://code.earthengine.google.com/628832ffbd371d828a048e9de2f0903d)). The energy balance output variables (`Hc`, `Hs`, `LEc`, `LEs`, `G`, `Rns`, and `Rnc`) obtained using GEE were also exported for comparison. \n",
    "\n",
    "Here we run the model and compare the output from GEE and our local run using xarray. \n",
    "\n",
    "## Reading the data\n",
    "\n",
    "Here we can read the data directly from the url, or you can download it and then replace the url with the local file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rioxarray\n",
    "\n",
    "url = \"https://raw.githubusercontent.com/kaust-halo/geeet/main/examples/data/landsat_era5_sample.tif\"\n",
    "data = rioxarray.open_rasterio(url, band_as_variable=True)\n",
    "data = (data.rename({band:data[band].attrs[\"long_name\"] for band in data})\n",
    "    .assign(\n",
    "        LE = lambda x: x.LEc+x.LEs,\n",
    "        H = lambda x: x.Hc+x.Hs,\n",
    "        Rn = lambda x: x.Rnc+x.Rns \n",
    "    )\n",
    ")\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's keep only the inputs for our local model run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs = data[[\"NDVI\", \"albedo\", \"radiometric_temperature\",  \n",
    "\"surface_pressure\", \"air_temperature\", \"dewpoint_temperature\", \"wind_speed\", \"solar_radiation\", \"thermal_radiation\"]]\n",
    "inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GEE ET model uses [ee.Image.pixelLonLat](https://developers.google.com/earth-engine/apidocs/ee-image-pixellonlat) internally whenever longitude and latitude information are needed. For the local ET model run, `lon` and `lat` need to be provided. Currently, there are two options. Either we can use `.rio.reproject(\"EPSG:4326\")` to transform the `x` and `y` coordinates to lon/lat, and then provide these, or we can provide a single `lon`, `lat` coordinate. For this small region, let's choose the latter. \n",
    "\n",
    "We also need to specify other scalar information including: doy (day of year), time, Vza (viewing zenith angle), zU (wind speed measurement height), and zT (temperature measurement height). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scalar_inputs = dict(\n",
    "    doy = 111,          # 2023-04-21; Day of year\n",
    "    time = 10.98,       # Local observation time\n",
    "    Vza = 0,            # Viewing zenith angle (degrees)\n",
    "    longitude = 38.25,  # Longitude (degrees)\n",
    "    latitude = 30.25,   # Latitude (degrees)\n",
    "    zU = 10,            # Wind measurement height (m)\n",
    "    zT = 2              # Temperature measurement height (m)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, run the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geeet\n",
    "\n",
    "xet = geeet.tseb.tseb_series(inputs, **scalar_inputs)\n",
    "xet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison\n",
    "\n",
    "First let's take an overview within the 0-350 W/m² range. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "xdiff = xet - data\n",
    "\n",
    "flux_name = \"LE\"\n",
    "plot_kwargs = dict(cmap=\"viridis\", vmin=0, vmax=350)\n",
    "\n",
    "facets = (xr.concat([\n",
    "    data[flux_name],\n",
    "    xet[flux_name],\n",
    "    xdiff[flux_name]\n",
    "], dim=\"band\")\n",
    ".assign_coords({\"band\": [\"GEE model\", \"xarray model\", \"difference\"]})\n",
    ".assign_attrs({\"units\": \"W/m²\"})\n",
    ")\n",
    "facetGrid = facets.plot(col=\"band\", **plot_kwargs)\n",
    "facetGrid.set_titles(template='{value}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further inspecting at a smaller LE scale, we can see the effect of using a single longitude/latitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "facets.sel(band=\"difference\").plot(vmin=-0.1, vmax=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reprojecting\n",
    "\n",
    "Let's try now using variable longitude and latitude. For this, we need to add them as variables to our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "datab = data.rio.reproject(\"EPSG:4326\")\n",
    "xv, yv = np.meshgrid(datab.x, datab.y, indexing='xy')\n",
    "lon, lat = np.meshgrid(datab.x, datab.y, indexing='xy')\n",
    "datab = datab.assign(longitude=([\"y\",\"x\"], lon), latitude=([\"y\",\"x\"], lat))\n",
    "\n",
    "inputsb = (datab[[\n",
    "    \"longitude\", \"latitude\",\n",
    "    \"NDVI\", \"albedo\", \"radiometric_temperature\", \n",
    "    \"surface_pressure\", \"air_temperature\", \"dewpoint_temperature\",\n",
    "    \"wind_speed\", \"solar_radiation\", \"thermal_radiation\"]]\n",
    ")\n",
    "\n",
    "xetb = geeet.tseb.tseb_series(inputsb, **dict(\n",
    "    doy = 111,          # 2023-04-21; Day of year\n",
    "    time = 10.98,       # Local observation time\n",
    "    Vza = 0,            # Viewing zenith angle (degrees)\n",
    "    zU = 10,            # Wind measurement height (m)\n",
    "    zT = 2              # Temperature measurement height (m)\n",
    "))\n",
    "\n",
    "xdiffb = xetb - datab # difference with the GEE model run\n",
    "\n",
    "xdiffb.LE.plot(vmin=-0.1, vmax=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the differences are no longer longitude-dependent, and much smaller in magnitude. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}