raster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0c2bc39d-4bf2-4bfc-8dc5-1978803d407a",
   "metadata": {},
   "source": [
    "# Remote Sensing Data\n",
    "\n",
    "This notebook works with some [Landsat 8 Collection 2 Level 2](https://planetarycomputer.microsoft.com/dataset/landsat-8-c2-l2) data, hosted by the Planetary Computer. We'll use the Planetary Computer's [STAC endpoint](https://planetarycomputer.microsoft.com/docs/reference/stac/) to query for a specific date range and area of interest, and [stackstac](https://stackstac.readthedocs.io/en/latest/) to construct an [xarray](http://xarray.pydata.org/) DataArray with all of our data.\n",
    "\n",
    "We'll use [Dask](https://dask.org/) to distribute the computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be3d4c49-ea8e-4828-ae44-d23b4bdd94e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import rioxarray\n",
    "\n",
    "import rasterio\n",
    "import stackstac\n",
    "import pystac_client\n",
    "import planetary_computer\n",
    "import rasterio.features\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs\n",
    "import matplotlib.animation\n",
    "import pandas as pd\n",
    "import xrspatial.multispectral\n",
    "from IPython.display import GeoJSON\n",
    "\n",
    "from dask_gateway import GatewayCluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5abda848-2e9f-4423-ba29-9c744850d2af",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = GatewayCluster()\n",
    "client = cluster.get_client()\n",
    "\n",
    "cluster.scale(72)\n",
    "cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37f3cac8-631a-4d79-b2b0-9b2b5741936b",
   "metadata": {},
   "outputs": [],
   "source": [
    "aoi = {\n",
    "  \"type\": \"FeatureCollection\",\n",
    "  \"features\": [\n",
    "    {\n",
    "      \"type\": \"Feature\",\n",
    "      \"properties\": {},\n",
    "      \"geometry\": {\n",
    "        \"type\": \"Polygon\",\n",
    "        \"coordinates\": [\n",
    "          [\n",
    "            [-93.0816650390625, 40.60144147645398],\n",
    "            [-91.724853515625, 40.60144147645398],\n",
    "            [-91.724853515625, 41.68111756290652],\n",
    "            [-93.0816650390625, 41.68111756290652],\n",
    "            [-93.0816650390625, 40.60144147645398]\n",
    "          ]\n",
    "        ]\n",
    "      }\n",
    "    }\n",
    "  ]\n",
    "}\n",
    "bbox = rasterio.features.bounds(aoi)\n",
    "GeoJSON(aoi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f638ba96-40a3-4bab-bd7d-48de956b6554",
   "metadata": {},
   "source": [
    "With `pystac_client` we can translate our desired \"give me this data covering this area over this time period\" into code quite naturally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4d043ca-fb2c-43af-8458-a7f720ccf080",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "stac = pystac_client.Client.open(\"https://planetarycomputer.microsoft.com/api/stac/v1\")\n",
    "\n",
    "search = stac.search(                 # I want...\n",
    "    collections=[\"landsat-8-c2-l2\"],  # this data...\n",
    "    bbox=bbox,                        # covering this area...\n",
    "    datetime=\"2016-01-01/2020-12-31\", # over this time period\n",
    "    limit=500,  # batch size\n",
    ")\n",
    "items = list(search.items())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd746f61-2381-40ec-9f9e-563fe1ce53e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "search.matched()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bedbf262-11f6-47ac-bb3f-2406bee6a624",
   "metadata": {},
   "outputs": [],
   "source": [
    "items[0].assets['SR_B3']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b1e1457-fe84-4db9-9087-d02fec252802",
   "metadata": {},
   "outputs": [],
   "source": [
    "items[0].properties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdb4d22f-e972-4185-b9d4-4c88f78208b0",
   "metadata": {},
   "source": [
    "The actual items in blob storage aren't publicly readable (we don't want people moving all of Landsat 8 out of Azure West Europe for \"free\"). You just need to \"sign\" the assets. You can sign them anonymously, you'll just be throttled. See https://planetarycomputer.microsoft.com/docs/concepts/sas/ for more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaa3b21e-32fb-44b2-a67c-cbebeea4ade1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "items = pystac_client.ItemCollection([\n",
    "    planetary_computer.sign_assets(item)\n",
    "    for item in items\n",
    "    if item.properties[\"eo:cloud_cover\"] < 25  # percent cloudy\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6dbc7e8-b689-4334-8041-ac3c7a4628b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ds = (\n",
    "    stackstac.stack(items, assets=[\"SR_B2\", \"SR_B3\", \"SR_B4\", \"SR_B5\"],\n",
    "                    bounds_latlon=bbox, chunksize=2056)\n",
    "        .assign_coords(band=[\"blue\", \"green\", \"red\", \"nir\"])\n",
    ").persist()\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e56aac9d-c6c9-41d8-b665-6150493dd662",
   "metadata": {},
   "outputs": [],
   "source": [
    "median = ds.median(dim=\"time\").persist()\n",
    "median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d5a571-9fde-4f90-8126-f18244706d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = xrspatial.multispectral.true_color(\n",
    "    *median.isel(x=slice(2000, 3000),\n",
    "                 y=slice(2000, 3000)).sel(band=[\"red\", \"green\", \"blue\"])\n",
    ").compute()\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 16),\n",
    "                       subplot_kw=dict(projection=cartopy.crs.AlbersEqualArea()))\n",
    "\n",
    "image = img.plot.imshow(ax=ax, add_labels=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3adc7dc1-05ac-42a3-8475-c2e519457088",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = xrspatial.ndvi(median.sel(band=\"nir\"), median.sel(band=\"red\")).persist()\n",
    "ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09b349cb-6e76-4c8d-8ce2-42bee8231774",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipyleaflet\n",
    "\n",
    "m = ipyleaflet.Map(center=[41.01447619471347, -92.48968322606777], scroll_wheel_zoom=True)\n",
    "m.layout.height = \"600px\"\n",
    "stackstac.add_to_map(ndvi, m, name=\"ndvi\", range=(-0.15, 0.5))\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20fff68f-9127-4f3e-8fa0-016e4d237f9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37e708e7-13aa-493c-ad4a-5fcb282834a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "red = ds.sel(band=\"red\")\n",
    "nir = ds.sel(band=\"nir\")\n",
    "\n",
    "ndvi_ts = (nir - red) / (nir + red)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c80289ed-5f7b-4f86-88dd-d7bdc24de134",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 6))\n",
    "\n",
    "ndvi_ts.median(dim=[\"y\", \"x\"]).plot(ax=ax);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}