05_EMIT_Waterbody_Extraction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# This is required while pygeos is installed\n",
    "os.environ['USE_PYGEOS'] = '0'\n",
    "import json\n",
    "\n",
    "import folium\n",
    "import geohash\n",
    "import geopandas\n",
    "import holoviews as hv\n",
    "import hvplot.xarray  # noqa: F401\n",
    "import numpy as np\n",
    "import odc.geo.xr\n",
    "import shapely\n",
    "import xarray as xr\n",
    "from dea_tools.spatial import xr_rasterize, xr_vectorize\n",
    "from fsspec.implementations.http import HTTPFileSystem\n",
    "from holoviews import opts\n",
    "\n",
    "from emit_tools import emit_xarray\n",
    "from utils import get_earthdata_token, get_rgb_dataset, load_emit_granule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Waterbody Extraction\n",
    "\n",
    "This notebook is a worked example. The goal of the analysis is to identify waterbodies that\n",
    "may have an algal bloom. To do this, we need to identify water bodies, examine their spectra\n",
    "and then calculate absorption depth before using a threshold to select waterbodies with high\n",
    "values.\n",
    "\n",
    "## Load data\n",
    "\n",
    "First we load data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# Loading data can take around 3-4 minutes on a 100 Mbps connection\n",
    "\n",
    "# See README.md for instructions on how to get an Earthdata token\n",
    "token = get_earthdata_token()\n",
    "\n",
    "# Refer to the README.md for instructions on how to find different granule IDs\n",
    "granule = \"EMIT_L2A_RFL_001_20230131T221923_2303114_008\"  # Kerang Lakes\n",
    "\n",
    "ds = load_emit_granule(granule, token)\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cleaning data\n",
    "\n",
    "Next we clean up the empty bands, before creating a simple water/not-water\n",
    "layer using the normalised-difference wetness index and bands that had\n",
    "high or low values over water as identified in the previous animations notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Clean up empty bands.\n",
    "ds = ds.fillna(np.nan).where(ds.reflectance!=-0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create a water layer\n",
    "high = ds.reflectance.sel(bands=450, method=\"nearest\")\n",
    "low = ds.reflectance.sel(bands=1275, method=\"nearest\")\n",
    "\n",
    "water = ((high - low) / (high + low)) > 0\n",
    "ds[\"water\"] = water.fillna(float(\"nan\")).where(water)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds.water.hvplot(aspect=\"equal\", frame_width=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering and creating IDs\n",
    "\n",
    "First, we want to filter by area so we only have large waterbodies. We\n",
    "do this using projected coordinates, rather than latitude/longitudes.\n",
    "\n",
    "Next we compute a geohash of the centroid of each waterbody, which can be\n",
    "used to help give a label to the waterbody, without needing to know the\n",
    "name of it.\n",
    "\n",
    "Note that we shrink the waterbody by a small amount."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "MIN_AREA = 80  # Hectares\n",
    "\n",
    "\n",
    "def add_geohash(row):\n",
    "    return geohash.encode(row.geometry.centroid.y, row.geometry.centroid.x, precision=9)\n",
    "\n",
    "\n",
    "# Create polygons from the water layer\n",
    "water_polygons = xr_vectorize(ds.water, crs=\"epsg:4326\", mask=ds.water.values == 1)\n",
    "water_polygons[\"area\"] = water_polygons.to_crs(\"epsg:3577\").area / 10000\n",
    "\n",
    "# Drop geopandas rows where the area is less than MIN_AREA\n",
    "water_polygons = water_polygons.drop(\n",
    "    water_polygons[water_polygons[\"area\"] < MIN_AREA].index\n",
    ")\n",
    "\n",
    "# Compute a geohash for each polygon at level 9\n",
    "geohashes = []\n",
    "for _, row in water_polygons.iterrows():\n",
    "    geohashes.append(add_geohash(row))\n",
    "\n",
    "water_polygons[\"geohash\"] = geohashes\n",
    "\n",
    "# Add an ID row\n",
    "water_polygons[\"id\"] = range(1, water_polygons.shape[0] + 1)\n",
    "\n",
    "# Show us what we've got\n",
    "print(\n",
    "    f\"Found {water_polygons.shape[0]} water polygons that are larger than {MIN_AREA} hectare(s)\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# View the water layer on an interactive map\n",
    "\n",
    "# Reduce the polygons by a small amount (in meters)\n",
    "SHRINK_AMOUNT = 100\n",
    "\n",
    "m = folium.Map(control_scale=True, tiles=None)\n",
    "\n",
    "for _, row in water_polygons.iterrows():\n",
    "    geometry = water_polygons.geometry.to_crs(\"epsg:3577\").buffer(-1 * SHRINK_AMOUNT).to_crs(\"epsg:4326\")\n",
    "    geojson = folium.GeoJson(\n",
    "        data=json.dumps(\n",
    "            shapely.geometry.mapping(geometry)\n",
    "        ),\n",
    "        style_function=lambda x: {\"fillColor\": \"blue\", \"Color\": \"blue\"},\n",
    "        tooltip=f\"{row.geohash}\",\n",
    "    )\n",
    "    folium.Popup(\n",
    "        f\"<p><strong>geohash:</strong> {row.geohash}<br><strong>area:</strong> {row['area']:.3f} Ha</p>\"\n",
    "    ).add_to(geojson)\n",
    "    geojson.add_to(m)\n",
    "\n",
    "# Zoom map\n",
    "m.fit_bounds(ds.odc.map_bounds())\n",
    "\n",
    "tile = folium.TileLayer(\n",
    "    tiles=\"https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}\",\n",
    "    attr=\"Esri\",\n",
    "    name=\"Esri Satellite\",\n",
    "    control=True,\n",
    ").add_to(m)\n",
    "\n",
    "folium.LayerControl().add_to(m)\n",
    "display(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Rasterise the polygons again, so we can join on the geohash later\n",
    "# First shrink them by SHRINK_AMOUNT\n",
    "water_polygons.geometry = water_polygons.geometry.to_crs(\"epsg:3577\").buffer(-1 * SHRINK_AMOUNT).to_crs(\"epsg:4326\")\n",
    "water_raster = xr_rasterize(water_polygons, ds, attribute_col=\"id\", crs=\"epsg:4326\")\n",
    "\n",
    "# Join the rasterised polygons to the dataset\n",
    "ds[\"id\"] = xr.DataArray(water_raster, dims=(\"latitude\", \"longitude\"))\n",
    "\n",
    "# Create another empty array of strings\n",
    "ds[\"geohash\"] = xr.DataArray(\n",
    "    np.full((ds.latitude.size, ds.longitude.size), \"\", dtype=\"U9\"),\n",
    "    dims=(\"latitude\", \"longitude\"),\n",
    ")\n",
    "\n",
    "for _, row in water_polygons.iterrows():\n",
    "    # I think 'where' works the opposite of what you'd expect\n",
    "    ds[\"geohash\"] = ds.geohash.where(ds.id != row.id, row.geohash)\n",
    "\n",
    "# Mask the empty values\n",
    "ds[\"geohash\"] = ds.geohash.where(ds.geohash != \"\", drop=False)\n",
    "del ds[\"id\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate summaries\n",
    "\n",
    "This next section creates summary values per waterbody (geohash).\n",
    "\n",
    "These will help us understand the average response and variability of\n",
    "the waterbody."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%capture --no-stdout\n",
    "\n",
    "means = ds.groupby(\"geohash\").mean()\n",
    "std_dev = ds.groupby(\"geohash\").std()\n",
    "min = ds.groupby(\"geohash\").min()\n",
    "max = ds.groupby(\"geohash\").max()\n",
    "\n",
    "# Create a new dataset with the mean, standard deviation, min and max values\n",
    "# for each geohash\n",
    "water_summaries = xr.Dataset(\n",
    "    {\n",
    "        \"mean\": means.reflectance,\n",
    "        \"std_dev\": std_dev.reflectance,\n",
    "        \"min\": min.reflectance,\n",
    "        \"max\": max.reflectance,\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Mean and std_dev plots\n",
    "color_cycle = hv.Cycle(\"Category20\")\n",
    "\n",
    "plots = []\n",
    "for geohash in water_summaries.geohash.values:\n",
    "    row = water_summaries.sel(geohash=geohash)\n",
    "\n",
    "    plots.append(\n",
    "        (\n",
    "            hv.Spread(row, vdims=[\"mean\", \"std_dev\", \"std_dev\"], label=f\"{geohash}\")\n",
    "            * hv.Curve(row, vdims=\"mean\", label=f\"{geohash}\")\n",
    "        )\n",
    "    )\n",
    "\n",
    "hv.Layout(plots).opts(\n",
    "    opts.Spread(color=color_cycle, show_legend=True),\n",
    "    opts.Curve(color=color_cycle, show_legend=True),\n",
    "    opts.Overlay(\n",
    "        show_title=True, frame_width=200, frame_height=50, show_legend=False, yaxis=None\n",
    "    ),\n",
    ").cols(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalised Difference Chlorophyll Index\n",
    "\n",
    "Reference: [https://doi.org/10.1016/j.rse.2011.10.016](https://doi.org/10.1016/j.rse.2011.10.016).\n",
    "\n",
    "The Normalised Difference Chlorophyll Index can be calculated using either Sentinel-2 or EMIT data\n",
    "Now we calculate the NDCI for each pixel for each waterbody.  This means we can use violin plots to further\n",
    "understand the variability of NDCI within each waterbody."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Calculate NDCI\n",
    "reference_band = 700\n",
    "absorption_band = 671\n",
    "\n",
    "ds[\"reference_band\"] = ds.reflectance.sel(bands=reference_band, method = 'nearest')\n",
    "ds[\"absorption_band\"] = ds.reflectance.sel(bands=absorption_band, method = 'nearest')\n",
    "ds[\"ND\"] = ((ds.reference_band - ds.absorption_band)/(ds.reference_band + ds.absorption_band))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot, masking to show only the water areas we identified earlier\n",
    "ds.ND.where(ds.water == 1).hvplot(aspect=\"equal\", frame_width=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Simplify to a summary dataset, removing the bands dimension and reflectance data\n",
    "ds_summary = ds.drop_dims(\"bands\")\n",
    "\n",
    "# Violin plots grouped by geohash\n",
    "ds_summary.hvplot.violin(\n",
    "    y=\"ND\",\n",
    "    by=\"geohash\",\n",
    ").opts(\n",
    "    opts.Violin(\n",
    "        width=1000,\n",
    "        height=600,\n",
    "        xrotation=45,\n",
    "        show_legend=False,\n",
    "        title=\"Normalised Difference Chlorophyll Index\",\n",
    "        ylim=(-0.1, 0.5),\n",
    "        violin_fill_color='ND',\n",
    "        cmap = 'Spectral_r',\n",
    "        clim = (0, 0.4),\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Absorption depth\n",
    "\n",
    "Reference: [https://doi.org/10.1080/01431161003789549](https://doi.org/10.1080/01431161003789549).\n",
    "\n",
    "Now we calculate the absorption depth for each waterbody. This is done\n",
    "on the underlying reflectance values over three different bands, and not\n",
    "on the summary data. This means we can use violin plots to further\n",
    "understand the variability of this new index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Calculate absorption depth, to be plotted per geohash\n",
    "absorption = 627\n",
    "reference_band = 560\n",
    "reference_band2 = 648\n",
    "\n",
    "absorption = ds.reflectance.sel(bands=absorption, method = 'nearest')\n",
    "reference1= ds.reflectance.sel(bands=reference_band, method = 'nearest')\n",
    "reference2 = ds.reflectance.sel(bands=reference_band2, method = 'nearest')\n",
    "ds[\"absorption_depth\"] = (reference1 + reference2)/2 - absorption"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Simplify to a summary dataset, removing the bands dimension and reflectance data\n",
    "ds_summary_2 = ds.drop_dims(\"bands\")\n",
    "\n",
    "# Violin plots grouped by geohash\n",
    "ds_summary_2.hvplot.violin(\n",
    "    y=\"absorption_depth\",\n",
    "    by=\"geohash\",\n",
    ").opts(\n",
    "    opts.Violin(\n",
    "        width=1000,\n",
    "        height=600,\n",
    "        xrotation=45,\n",
    "        show_legend=False,\n",
    "        title=\"Absorption Depth\",\n",
    "        ylim=(-0.02, 0.03),\n",
    "        violin_fill_color='absorption_depth',\n",
    "        cmap = 'Spectral_r',\n",
    "        clim = (0, 0.02),\n",
    "    )\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Picking a threshold\n",
    "\n",
    "Here we're going to select water bodies (geohashes) by thresholding\n",
    "the absorption depth. First plot the absorption depths on a scatter plot.\n",
    "The x axis is just the waterbody's geohash, so it's meaningless.\n",
    "\n",
    "The following cell uses a threshold of 0.01, so selects the three\n",
    "waterbodies with a mean absorption_depth of above this and plots the\n",
    "spectra for these three waterbodies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "geohash_absorption_depths = ds_summary_2.groupby(\"geohash\").mean()\n",
    "\n",
    "geohash_absorption_depths.hvplot(\n",
    "    x=\"geohash\",\n",
    "    y=\"absorption_depth\",\n",
    "    kind=\"scatter\",\n",
    "    title=\"Mean absorption depth\",\n",
    "    color=\"absorption_depth\",\n",
    "    cmap = 'magma_r',\n",
    "    # Line color\n",
    "    line_color=\"grey\",\n",
    "    size=40,\n",
    "    xaxis=None,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "abs_d_gt_001 = geohash_absorption_depths.where(\n",
    "    geohash_absorption_depths.absorption_depth > 0.01, drop=True\n",
    ")\n",
    "high_absv = list(abs_d_gt_001.geohash.values)\n",
    "\n",
    "# Mean and std_dev plots\n",
    "color_cycle = hv.Cycle(\"Category20\")\n",
    "\n",
    "plots = []\n",
    "for geohash in high_absv:\n",
    "    row = water_summaries.sel(geohash=geohash)\n",
    "\n",
    "    plots.append(\n",
    "        (\n",
    "            hv.Spread(row, vdims=[\"mean\", \"std_dev\", \"std_dev\"], label=f\"{geohash}\")\n",
    "            * hv.Curve(row, vdims=\"mean\", label=f\"{geohash}\")\n",
    "        )\n",
    "    )\n",
    "\n",
    "# Add a mean of all waterbodies plot too\n",
    "mean_all = water_summaries.mean(\"geohash\")\n",
    "plots.append(\n",
    "    (\n",
    "        hv.Spread(mean_all, vdims=[\"mean\", \"std_dev\", \"std_dev\"], label=f\"all\")\n",
    "        * hv.Curve(mean_all, vdims=\"mean\", label=f\"all\")\n",
    "    )\n",
    ")\n",
    "\n",
    "hv.Overlay(plots).opts(\n",
    "    opts.Spread(color=color_cycle, show_legend=True),\n",
    "    opts.Curve(color=color_cycle, show_legend=True),\n",
    "    opts.Overlay(\n",
    "        show_title=True, frame_width=600, frame_height=300, show_legend=True, yaxis=None\n",
    "    ),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show the selected waterbodies on a map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m = folium.Map(control_scale=True, tiles=None)\n",
    "\n",
    "for _, row in water_polygons.iterrows():\n",
    "    # Skip the ones we haven't selected\n",
    "    if row.geohash not in high_absv:\n",
    "        continue\n",
    "    geojson = folium.GeoJson(\n",
    "        data=json.dumps(\n",
    "            shapely.geometry.mapping(row.geometry)\n",
    "        ),\n",
    "        style_function=lambda x: {\"fillColor\": \"blue\", \"Color\": \"blue\"},\n",
    "        tooltip=f\"{row.geohash}\",\n",
    "    )\n",
    "    folium.Popup(\n",
    "        f\"<p><strong>geohash:</strong> {row.geohash}<br><strong>area:</strong> {row['area']:.3f} Ha</p>\"\n",
    "    ).add_to(geojson)\n",
    "    geojson.add_to(m)\n",
    "\n",
    "# Zoom map\n",
    "m.fit_bounds(ds.odc.map_bounds())\n",
    "\n",
    "tile = folium.TileLayer(\n",
    "    tiles=\"https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}\",\n",
    "    attr=\"Esri\",\n",
    "    name=\"Esri Satellite\",\n",
    "    control=True,\n",
    ").add_to(m)\n",
    "\n",
    "folium.LayerControl().add_to(m)\n",
    "display(m)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Alex Test Environment",
   "language": "python",
   "name": "test"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}