notebooks/odc-landsat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5071afec",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Landsat Collection 2 on AWS\n",
    "\n",
    "This notebook explores the Landsat Collection 2 data on AWS:\n",
    "\n",
    "- [Landsat Collection 2 STAC API](https://landsatlook.usgs.gov/stac-server), a catalog of Landsat data\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",
    "- [Coiled.io](https://coiled.io/), a service for hosting Dask clusters\n",
    "- [hvplot](https://hvplot.holoviz.org/) for visualization\n",
    "\n",
    "Shown will be how to 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": "3120189b-b49c-41d5-a702-88548ba24b49",
   "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",
    "import os\n",
    "os.environ['AWS_REQUEST_PAYER'] = 'requester'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac9e6132-9bc6-4a97-a00d-b2f40159251d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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": "0f100b1b",
   "metadata": {},
   "source": [
    "# Search for data\n",
    "\n",
    "Use pystac-client to find data in the Landsat STAC API. First, print up a table of all the STAC Collections in the API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d64615c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the Landsat STAC API\n",
    "\n",
    "from pystac_client import Client\n",
    "URL = 'https://landsatlook.usgs.gov/stac-server'\n",
    "cat = Client.open(URL)\n",
    "print(cat)\n",
    "\n",
    "collections = [(c.id, c.title) for c in cat.get_collections()]\n",
    "pd.set_option(\"display.max_colwidth\", 150)\n",
    "df = pd.DataFrame(collections, columns=['id', 'title'])\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7b0b19a",
   "metadata": {},
   "source": [
    "Fetch the collection of interest: Landsat Collection 2, Level 2 Surface Reflectance (landsat-c2l2-sr) and print the assets that are available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21053106",
   "metadata": {},
   "outputs": [],
   "source": [
    "collection_id = 'landsat-c2l2-sr'\n",
    "\n",
    "collection = cat.get_collection(collection_id)\n",
    "pd.DataFrame.from_dict(collection.to_dict()['item_assets'], orient='index')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "972e0654",
   "metadata": {},
   "source": [
    "Change the AOI, search parameters here, and print how many matching scenes there are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42062ad1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import json\n",
    "\n",
    "aoi = gpd.read_file('../aois/bear-fire.geojson')\n",
    "geom = json.loads(aoi['geometry'].to_json())['features'][0]['geometry']\n",
    "\n",
    "# limit sets the # of items per page so we can see multiple pages getting fetched\n",
    "search = cat.search(\n",
    "    collections = [collection_id],\n",
    "    intersects = aoi['geometry'][0],\n",
    "    datetime = \"2021-01-01/2021-03-31\",\n",
    "    query = [\"eo:cloud_cover<10\"],\n",
    "    limit = 100\n",
    ")\n",
    "\n",
    "print(f\"{search.matched()} items found\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "734177ef",
   "metadata": {},
   "source": [
    "# Use GeoPandas to view footprints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "946997a9",
   "metadata": {},
   "source": [
    "The cell below fetches all the STAC Items and updates the URLs to use the provided s3 URLs which can be used for direct access rather than the default https URLs. This is because the alternate extension is not yet supported in PySTAC, when it is there will be an easier way to specify which alternate URL, if any, to use for the assets.\n",
    "\n",
    "Then, we create a GeoDataFrame for visualizing the footprints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fa5182f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all items as a dictionary\n",
    "items_dict = search.get_all_items_as_dict()['features']\n",
    "\n",
    "# update URLs to use s3\n",
    "for item in items_dict:\n",
    "    for a in item['assets']:\n",
    "        if 'alternate' in item['assets'][a] and 's3' in item['assets'][a]['alternate']:\n",
    "            item['assets'][a]['href'] = item['assets'][a]['alternate']['s3']['href']\n",
    "        item['assets'][a]['href'] = item['assets'][a]['href'].replace('usgs-landsat-ard', 'usgs-landsat')\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",
    "pd.reset_option(\"display.max_colwidth\")\n",
    "items_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e37d05f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_polygons(aoi) * items_gdf.hvplot.paths(geo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f258702",
   "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": "74858c15",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "\n",
    "cfg = \"\"\"---\n",
    "landat-c2l2-sr:\n",
    "  measurements:\n",
    "    '*':\n",
    "      dtype: uint16\n",
    "      nodata: 0\n",
    "      unit: 'm'\n",
    "\"\"\"\n",
    "cfg = yaml.load(cfg, Loader=yaml.CSafeLoader)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08311f41",
   "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, chunk size, and geometry to clip the data to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "583cd3a8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "from odc.stac import stac_load\n",
    "\n",
    "# Create PySTAC ItemCollection\n",
    "item_collection = pystac.ItemCollection(items_dict)\n",
    "\n",
    "# default to CRS and resolution from first Item\n",
    "from pystac.extensions.projection import ProjectionExtension\n",
    "from pyproj import CRS\n",
    "\n",
    "dc = stac_load(item_collection,\n",
    "               bands=['red', 'blue', 'green', 'nir08'],\n",
    "               chunks={\"x\": 2048, \"y\": 2048},\n",
    "               groupby='solar_day',\n",
    "               stac_cfg=cfg,\n",
    "               geopolygon=geom,\n",
    ")\n",
    "dc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8920daa1",
   "metadata": {},
   "source": [
    "# Calculations\n",
    "\n",
    "We create an RGB datacube representation and generate an NDVI datacube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0adc4f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "vis = dc.odc.to_rgba(vmin=1, vmax=20000, bands=['red', 'green', 'blue'])\n",
    "vis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ed4192d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = ((dc['nir08'] - dc['red']) / (dc['nir08'] + dc['red'])).clip(0, 1)\n",
    "ndvi.name = 'ndvi'\n",
    "ndvi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99b11cdb",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Start Dask Client\n",
    "\n",
    "Start either a local Dask, or use [coiled.io](coiled.io)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b51769f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# local Dask\n",
    "\n",
    "from dask.distributed import Client\n",
    "client = Client()\n",
    "client"
   ]
  },
  {
   "cell_type": "raw",
   "id": "23d1937b",
   "metadata": {},
   "source": [
    "%%time\n",
    "\n",
    "import coiled\n",
    "from dask.distributed import Client\n",
    "\n",
    "# start dask cluster on coiled.io\n",
    "cluster = coiled.Cluster(\n",
    "    n_workers=25, software=\"matthewhanson/geo-notebooks-basic\", backend_options={\"region\": \"us-west-2\"},\n",
    "    environ={\"GDAL_DISABLE_READDIR_ON_OPEN\": \"YES\", \"AWS_REQUEST_PAYER\": \"requester\"}\n",
    ")\n",
    "# re-use existing cluster\n",
    "#cluster = coiled.Cluster(name=\"matthewhanson-e6818fa6-e\")\n",
    "client = Client(cluster)\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b22953d",
   "metadata": {
    "tags": []
   },
   "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": "c45736e3-9d25-4660-b29d-816114e88472",
   "metadata": {},
   "outputs": [],
   "source": [
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98c54c0c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "from dask.distributed import wait\n",
    "\n",
    "vis = client.persist(vis)\n",
    "_ = wait(vis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae6b5594",
   "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": "21d06098-4a08-48f0-beb0-5c83f19f6aae",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.xarray\n",
    "import panel as pn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f400ecef-fd6f-4126-89b3-802a5bf57fc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "vis_.spatial_ref.crs_wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acfd80ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "hvplot_kwargs = dict(crs=32610, tiles='OSM', \n",
    "                frame_width=400, widgets={'time': pn.widgets.Select})\n",
    "\n",
    "vis_.hvplot.rgb(x='x', y='y', bands='band', **hvplot_kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f90dcb0-bcfb-4604-adde-92f34a18ce5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "vis_.hvplot.rgb(x='x', y='y', bands='band', **hvplot_kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ad16c61",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ndvi_ = ndvi.compute()\n",
    "ndvi_.hvplot(x='x', y='y', **hvplot_kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6326865",
   "metadata": {},
   "source": [
    "Create an animated GIF of NDVI over time using `geogif` with the fetched results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "028dafda",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geogif import gif\n",
    "\n",
    "#ndvi_ = ndvi_c.transpose('time','x','y').compute()\n",
    "\n",
    "gif(ndvi_, fps=1, cmap='YlGn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6c27e1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ndvi_mean = ndvi.mean(dim=['x', 'y']).compute()\n",
    "ndvi_mean.hvplot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb97978a",
   "metadata": {},
   "source": [
    "# Shutdown cluster\n",
    "\n",
    "Shut down the cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f21e268d",
   "metadata": {},
   "outputs": [],
   "source": [
    "client.close()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "odc:Python",
   "language": "python",
   "name": "conda-env-odc-py"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}