notebooks/odc-planetary-computer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "23fedf2a",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Sentinel-2 on the Planetary Computer\n",
    "\n",
    "This notebook explores Sentinel-2 data on Microsoft's Planetary Computer using:\n",
    "\n",
    "- [Planetary Computer STAC API](https://planetarycomputer.microsoft.com/api/stac/v1), catalog of public data\n",
    "- [Planetary Computer Hub](https://planetarycomputer.microsoft.com/) for running Jupyter Notebooks in the cloud\n",
    "- [pystac-client](https://pystac-client.readthedocs.io/) for searching and access data\n",
    "- [OpenDataCube](https://www.opendatacube.org/) and [odc-stac](https://odc-stac.readthedocs.io/) for loading STAC assets and representing geospatial data as XArrays\n",
    "- [XArray](http://xarray.pydata.org/en/stable/), [pandas](https://pandas.pydata.org/) and [geopandas](https://geopandas.org/) for manipulating data\n",
    "- [Dask](https://dask.org/) for performing parallel, distributed computing\n",
    "- [hvplot](https://hvplot.holoviz.org/) for visualization\n",
    "\n",
    "Shown will be how find data for an area of interest, explore the resulting metadata, perform calculations, and visualize the results.\n",
    "\n",
    "Created by [Element 84](http://element84.com/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ee670d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial imports and reusable functions\n",
    "\n",
    "import holoviews as hv\n",
    "hv.extension('bokeh')\n",
    "\n",
    "from copy import deepcopy\n",
    "import geopandas as gpd\n",
    "import hvplot.pandas\n",
    "import pandas as pd\n",
    "import pystac\n",
    "from shapely.geometry import shape\n",
    "\n",
    "# create a function for later reuse\n",
    "def plot_polygons(data, *args, **kwargs):\n",
    "    return data.hvplot.paths(*args, geo=True, tiles='OSM', xaxis=None, yaxis=None,\n",
    "                             frame_width=600, frame_height=600,\n",
    "                             line_width=3, **kwargs)\n",
    "\n",
    "# convert a list of STAC Items into a GeoDataFrame\n",
    "def items_to_geodataframe(items):\n",
    "    _items = []\n",
    "    for i in items:\n",
    "        _i = deepcopy(i)\n",
    "        _i['geometry'] = shape(_i['geometry'])\n",
    "        _items.append(_i)\n",
    "    gdf = gpd.GeoDataFrame(pd.json_normalize(_items))\n",
    "    for field in ['properties.datetime', 'properties.created', 'properties.updated']:\n",
    "        if field in gdf:\n",
    "            gdf[field] = pd.to_datetime(gdf[field])\n",
    "    gdf.set_index('properties.datetime', inplace=True)\n",
    "    return gdf\n",
    "\n",
    "# set pystac_client logger to DEBUG to see API calls\n",
    "import logging\n",
    "logging.basicConfig()\n",
    "logger = logging.getLogger('pystac_client')\n",
    "logger.setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e5a8e17",
   "metadata": {},
   "source": [
    "# Search for data\n",
    "\n",
    "Use pystac-client to find data in the Planetary Computer STAC API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f67b8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the Planetary Computer STAC API\n",
    "\n",
    "from pystac_client import Client\n",
    "URL = 'https://planetarycomputer.microsoft.com/api/stac/v1'\n",
    "cat = Client.open(URL)\n",
    "cat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c204a110",
   "metadata": {},
   "source": [
    "Fetch the collection of interest: Sentinel-2, Level 2 Surface Reflectance and print the assets that are available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16b7e2c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "collection = cat.get_collection('sentinel-2-l2a')\n",
    "\n",
    "pd.DataFrame.from_dict(collection.to_dict()['item_assets'], orient='index')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8c1f8a6",
   "metadata": {},
   "source": [
    "Change the AOI, search parameters here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f187e5d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import json\n",
    "\n",
    "aoi = gpd.read_file('../aois/bear-fire.geojson')\n",
    "geom = aoi['geometry'][0] # < shapely geometry object\n",
    "\n",
    "# limit sets the # of items per page so we can see multiple pages getting fetched\n",
    "search = cat.search(\n",
    "    collections = [\"sentinel-2-l2a\"],\n",
    "    intersects = geom,\n",
    "    datetime = \"2019-10-01/2019-10-31\",\n",
    "    query = [\"eo:cloud_cover<25\"],\n",
    "    limit = 100\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1da0db1f",
   "metadata": {},
   "source": [
    "# Use GeoPandas to view footprints\n",
    "\n",
    "The cell below fetches all the STAC Items, then creates a GeoDataFrame for visualizing the footprints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6568c07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all items as a dictionary\n",
    "items_dict = search.get_all_items_as_dict()['features']\n",
    "\n",
    "# Create GeoDataFrame from Items\n",
    "items_gdf = items_to_geodataframe(items_dict)\n",
    "\n",
    "print(f\"{len(items_dict)} items found\")\n",
    "\n",
    "items_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1c9e06c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_polygons(aoi) * items_gdf.hvplot.paths(geo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe1dc063",
   "metadata": {},
   "source": [
    "# OpenDataCube\n",
    "\n",
    "Now we'll turn the set of scenes into a virtual datacube. None of the data will actually be read yet.\n",
    "\n",
    "The configuration string (`cfg`) is for providing additional info not currently available in the STAC Items, but will be in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb7e5483",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "\n",
    "cfg = \"\"\"---\n",
    "sentinel-2-l2a:\n",
    "  assets:\n",
    "    '*':\n",
    "      data_type: uint16\n",
    "      nodata: 0\n",
    "      unit: '1'\n",
    "\"\"\"\n",
    "cfg = yaml.load(cfg, Loader=yaml.CSafeLoader)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b912655",
   "metadata": {},
   "source": [
    "Here we load as a DataCube. A PySTAC ItemCollection is created from the found STAC Items, and we specify various parameters, such as bands of interest and chunk size. We are requesting to only load pixels within a bounding box of the requested geometry (`bbox=geom.bounds`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd4e2099",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "from odc.stac import stac_load\n",
    "import planetary_computer as pc\n",
    "\n",
    "# Create PySTAC ItemCollection\n",
    "item_collection = pystac.ItemCollection(items_dict)\n",
    "\n",
    "dc = stac_load(item_collection,\n",
    "               measurements=['B02', 'B03', 'B04', 'B08'],\n",
    "               chunks={\"x\": 2048, \"y\": 2048},\n",
    "               bbox=geom.bounds,\n",
    "               stac_cfg=cfg,\n",
    "               patch_url=pc.sign\n",
    ")\n",
    "dc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "777ce4df",
   "metadata": {},
   "source": [
    "# Calculations\n",
    "\n",
    "We will create an RGBA datacube representation (`nodata` values have `alpha=0`), and generate an NDVI datacube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b1ce170",
   "metadata": {},
   "outputs": [],
   "source": [
    "from odc.algo import to_rgba\n",
    "\n",
    "RGB = ['B04', 'B03', 'B02']\n",
    "vis = to_rgba(dc, clamp=(1, 3000), bands=RGB)\n",
    "vis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8182986c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = ((dc['B08'] - dc['B04']) / (dc['B08'] + dc['B04'])).clip(0, 1).rename(\"ndvi\")\n",
    "ndvi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63da01f8",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Start Dask Client\n",
    "\n",
    "Start a local Dask cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbf52e91",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client\n",
    "client = Client()\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d4ca4ca",
   "metadata": {},
   "source": [
    "# Compute\n",
    "\n",
    "Now, we kick off our Dask computation by using the Dask persist function, which will keep the data in memory on the cluster for faster access later.\n",
    "\n",
    "The Dask `compute` function is used when we actually want the data, such as displaying it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "374855ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from dask.distributed import wait\n",
    "\n",
    "ndvi, vis = client.persist([ndvi, vis])\n",
    "_ = wait([ndvi, vis])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "668c90a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "vis_ = vis.compute()\n",
    "vis_.plot.imshow(col='time', rgb='band', col_wrap=5, robust=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85377d34",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.xarray\n",
    "\n",
    "hvplot_kwargs = {\n",
    "    \"frame_width\": 800,\n",
    "    \"xaxis\": None,\n",
    "    \"yaxis\": None,\n",
    "    \"widget_location\": \"bottom\",\n",
    "    \"aspect\": len(vis.x)/len(vis.y)\n",
    "}\n",
    "\n",
    "vis_.hvplot.rgb('x', 'y', bands='band', groupby='time', **hvplot_kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "354c8eda",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_ = ndvi.compute()\n",
    "ndvi_.hvplot('x', 'y', groupby='time', **hvplot_kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2709270",
   "metadata": {},
   "source": [
    "Create an animated GIF of NDVI over time using `geogif` with the fetched results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a347dea",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from geogif import gif, dgif\n",
    "\n",
    "gif(ndvi_, fps=5, cmap='YlGn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09a31a96",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ndvi_mean = ndvi.mean(dim=['x', 'y']).compute()\n",
    "ndvi_mean.hvplot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9128e55",
   "metadata": {},
   "source": [
    "# Shutdown cluster\n",
    "\n",
    "Shut down the cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b470887",
   "metadata": {},
   "outputs": [],
   "source": [
    "if \"cluster\" in locals():\n",
    "    cluster.close()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}