04_EMIT_Interactive_Polygons.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "import holoviews as hv\n",
    "import hvplot.xarray  # noqa: F401\n",
    "import numpy as np\n",
    "import rioxarray\n",
    "import xarray as xr\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,\n",
    "                   hv_stream_to_rio_geometries, load_emit_granule)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive polygon explorer\n",
    "\n",
    "This notebook allows you to draw polygons over an RGB image and get\n",
    "a summary of the spectra of all pixels in that polygon. It'll\n",
    "plot the mean and standard deviation of each polygon.\n",
    "\n",
    "## Load data"
   ]
  },
  {
   "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_20230316T045211_2307503_006\" # Canberra\n",
    "\n",
    "ds = load_emit_granule(granule, token)\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%capture --no-stdout\n",
    "# Create an RGB dataset for rendering\n",
    "\n",
    "from utils import get_rgb_dataset\n",
    "\n",
    "# Select wavelengths to be displayed\n",
    "r = 2100\n",
    "g = 950\n",
    "b = 490\n",
    "\n",
    "# Select brightness - range between 0-1, higher values 'brighten' the whole scene\n",
    "brightness = 0.4\n",
    "\n",
    "## End configuration area ##\n",
    "\n",
    "# Get data structured for rendering\n",
    "ds_rgb = get_rgb_dataset(ds, [r, g, b], brightness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean up empty bands.\n",
    "ds = ds.fillna(np.nan).where(ds.reflectance!=-0.01)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the plots\n",
    "\n",
    "This example has a different cell for the capture of areas and the rendering\n",
    "of the graphs, because the graphs take a long time to create.\n",
    "\n",
    "In this map that is output from the following cell, use the\n",
    "\"Polygon Draw Tool\" to draw up to five polygons, then run the next\n",
    "cell to see the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Limit the number of drawn polygons\n",
    "POLY_LIMIT = 5\n",
    "\n",
    "color_cycle = hv.Cycle('Category10')\n",
    "colors = [color_cycle[i] for i in range(5)]\n",
    "\n",
    "# RGB image/map\n",
    "map = ds_rgb.hvplot.rgb(\n",
    "    x=\"longitude\", y=\"latitude\", bands=\"bands\", aspect=\"equal\", frame_width=600\n",
    ")\n",
    "\n",
    "# Set up a holoviews points array to enable plotting of the clicked points\n",
    "xmid = ds.longitude.values[int(len(ds.longitude) / 2)]\n",
    "ymid = ds.latitude.values[int(len(ds.latitude) / 2)]\n",
    "polygons = hv.Polygons(\n",
    "    [],\n",
    "    kdims=[\"xs\", \"ys\"],\n",
    ")\n",
    "\n",
    "polygons_stream = hv.streams.PolyDraw(\n",
    "    # data=polygons.columns(),\n",
    "    source=polygons,\n",
    "    num_objects=POLY_LIMIT,\n",
    "    styles={'fill_color': color_cycle.values[0:POLY_LIMIT]}\n",
    ")\n",
    "\n",
    "# Plot the Map and Dynamic Map side by side\n",
    "(map * polygons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Build a spectral plot for each of the drawn polygons.\n",
    "# This takes a while.\n",
    "\n",
    "plots = []\n",
    "\n",
    "geometries = hv_stream_to_rio_geometries(polygons_stream.data)\n",
    "export = []\n",
    "\n",
    "for i, geometry in enumerate(geometries):\n",
    "    data = ds.reflectance.rio.clip(geometry, drop=False)\n",
    "    hv_data = hv.Dataset(data, kdims=[\"bands\", \"latitude\", \"longitude\"], vdims=[\"reflectance\"])\n",
    "    agg = hv_data.aggregate(\"bands\", np.nanmean, spreadfn=np.nanstd)\n",
    "    \n",
    "    data = [i, json.dumps(geometry)] + list(agg.data.reflectance.values) + list(agg.data.reflectance_nanstd.values)\n",
    "    export.append(data)                                                    \n",
    "\n",
    "    plots.append(\n",
    "        (hv.Spread(agg) * hv.Curve(agg, label=f\"{i}\"))\n",
    "    )\n",
    "\n",
    "hv.Overlay(plots).opts(\n",
    "    opts.Spread(color=color_cycle),\n",
    "    opts.Curve(color=color_cycle),\n",
    "    opts.Overlay(show_title=False, frame_width=600, show_legend=False)\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Output results to a CSV\n",
    "\n",
    "This CSV has a lot of columns, but there's no other way to export the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import csv \n",
    "\n",
    "data = polygons_stream.data\n",
    "wavelengths = ds.wavelengths.values\n",
    "\n",
    "rows_headings = [[\"id\", \"geom\"] + [f\"{i}_mean\" for i in wavelengths] + [f\"{i}_std\" for i in wavelengths]]\n",
    "\n",
    "rows = rows_headings + export\n",
    "\n",
    "with open('polygon_data.csv', 'w') as f:\n",
    "    writer = csv.writer(f, delimiter=\";\")\n",
    "    writer.writerows(rows)"
   ]
  }
 ],
 "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}