01_EMIT_Basic.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import hvplot.xarray  # noqa: F401\n",
    "import xarray as xr\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": [
    "# Starter Notebook\n",
    "\n",
    "This notebook demonstrates some basic functionality that is useful in working with the EMIT data\n",
    "from NASA.\n",
    "\n",
    "## Loading data\n",
    "\n",
    "To load data, there's a function provided by NASA that will restructure the NetCDF file,\n",
    "which we can stream from their cloud data store.\n",
    "\n",
    "To discover scenes, you can browse NASA's STAC catalog in the Radiant Earth [STAC Browser](https://radiantearth.github.io/stac-browser/#/external//cmr.earthdata.nasa.gov/cloudstac/LPCLOUD/collections/EMITL2ARFL.v001).\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_20230316T045211_2307503_006\" # Canberra\n",
    "\n",
    "ds = load_emit_granule(granule, token)\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Render data as RGB\n",
    "\n",
    "This section will select three bands to visualise on the red, green and blue\n",
    "channels of an image and will stretch the data to make it easier to view.\n",
    "\n",
    "This plot is very simple and is a useful visual aid, but it's no longer data, it's a visualisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 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",
    "# Get data structured for rendering\n",
    "ds_rgb = get_rgb_dataset(ds, [r, g, b], brightness)\n",
    "\n",
    "# Plot the RGB data\n",
    "ds_rgb.hvplot.rgb(\n",
    "    x=\"longitude\", y=\"latitude\", bands=\"bands\", aspect=\"equal\", frame_width=600\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run normalised band index\n",
    "\n",
    "Here we're going to visualise some data and convert it into information in the form\n",
    "of a true/false value for whether or not a pixel is water.\n",
    "\n",
    "The logic here is to pick a band where the reflectance of a pixel over water is\n",
    "high and where it's low, to then isolate those bands and finally do an NDWI calculation.\n",
    "\n",
    "From [this reference](https://bmcecol.biomedcentral.com/articles/10.1186/s12898-019-0233-0/tables/2),\n",
    "it's suggested to use something like `860` for the first value and\n",
    "`1240`, `1640` or `2130` for the high value (short/long wavelengths). By inspecting bands,\n",
    "I've picked `650` for the first band, as it has high values over water and `1275` as it\n",
    "has low values over water. This should lead to water having values approaching 1, and\n",
    "non-water with values approaching -1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# High meaning bright, or high reflectance\n",
    "high = 450\n",
    "low = 1275\n",
    "\n",
    "ds[\"water_high\"] = ds.reflectance.sel(bands=high, method=\"nearest\")\n",
    "ds[\"water_low\"] = ds.reflectance.sel(bands=low, method=\"nearest\")\n",
    "ds[\"water_diff\"] = ds.water_high - ds.water_low\n",
    "\n",
    "ds[\"ndwi\"] = ds.water_diff / (ds.water_high + ds.water_low)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_options = {\n",
    "    \"aspect\": \"equal\",\n",
    "    \"frame_width\": 300,\n",
    "    \"xaxis\": \"bare\",\n",
    "    \"yaxis\": \"bare\",\n",
    "}\n",
    "\n",
    "image_high = ds.water_high.hvplot(**plot_options, title=f\"high: {high}\", clim=(0, 0.4), cmap=\"viridis\")\n",
    "image_low = ds.water_low.hvplot(**plot_options, title=f\"low: {low}\", clim=(0, 0.4), cmap=\"viridis\")\n",
    "image_diff = ds.water_diff.hvplot(**plot_options, title=\"difference\", clim=(-1, 1), cmap=\"RdBu\")\n",
    "image_ndwi = ds.ndwi.hvplot(**plot_options, title=\"NDWI\", clim=(-1, 1), cmap=\"bwy_r\")\n",
    "\n",
    "\n",
    "(image_high + image_low + image_diff + image_ndwi).cols(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot a histogram of the NDWI values, so we can pick a threshold\n",
    "ds.ndwi.hvplot.hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create the new data variable\n",
    "ds[\"threshold_water\"] = ds.ndwi > 0.0\n",
    "# Mask out the zero values\n",
    "ds[\"threshold_water\"] = ds.threshold_water.where(ds.threshold_water!=0)\n",
    "# Plot the result\n",
    "ds.threshold_water.hvplot(frame_width=600, aspect=\"equal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save and load xarray\n",
    "\n",
    "Saving work to disk as a cache or as an input to a separate process is simple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Store xarray dataset on disk. The Canberra example was 5 GB when saved.\n",
    "ds.to_netcdf(\"saved_on_disk.nc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load the xarray from disk again\n",
    "ds_new = xr.open_dataset(\"saved_on_disk.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exporting raster data\n",
    "\n",
    "Finally, we'll export data in a standard geotiff format, which can then be used in\n",
    "a GIS application or in web mapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Convert threshold_water band to uint8 and save to disk\n",
    "ds.threshold_water.astype(\"uint8\").rio.to_raster(f\"water.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Export a single band as a GeoTIFF\n",
    "ds.sel(bands=r, method=\"nearest\").reflectance.rio.to_raster(f\"band_{r}.tif\")\n",
    "ds.sel(bands=g, method=\"nearest\").reflectance.rio.to_raster(f\"band_{g}.tif\")\n",
    "ds.sel(bands=b, method=\"nearest\").reflectance.rio.to_raster(f\"band_{b}.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Combine three images into a single GeoTIFF with three bands\n",
    "from osgeo_utils import gdal_merge\n",
    "\n",
    "output_file_path = f\"bands_{r}_{g}_{b}.tif\"\n",
    "input_files_path = f\"band_{r}.tif band_{g}.tif band_{b}.tif\".split()\n",
    "parameters = ['', '-o', output_file_path] + input_files_path + ['-separate', '-co', 'COMPRESS=LZW']\n",
    "gdal_merge.main(parameters)"
   ]
  }
 ],
 "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}