datasets/chesapeake/chesapeake-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d2583da9-3f1e-40e6-b81d-700f26e89daa",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Accessing Chesapeake Land Use data with the Planetary Computer STAC API\n",
    "\n",
    "The [Chesapeake Land Use](https://www.chesapeakeconservancy.org/conservation-innovation-center/high-resolution-data/land-use-data-project/) dataset is high-resolution 1-meter Land Use data product in raster format for the entire Chesapeake Bay watershed. The dataset was created by modifying the 2013-2014 high-resolution land cover dataset using 13 ancillary datasets including data on zoning, land use, parcel boundaries, landfills, floodplains, and wetlands. The product area encompasses over 250,000 square kilometers in New York, Pennsylvania, Maryland, Delaware, West Virginia, Virginia, and the District of Columbia. The dataset was created by the Chesapeake Conservancy Conservation Innovation Center for the Chesapeake Bay Program, which is a regional partnership of EPA, other federal, state, and local agencies and governments, nonprofits, and academic institutions that leads and directs Chesapeake Bay restoration efforts.\n",
    "\n",
    "In this notebook, we'll demonstrate how to access and work with this data through the Planetary Computer. Documentation for this dataset is available at the [Planetary Computer Data Catalog](https://planetarycomputer.microsoft.com/dataset/chesapeake-lu).\n",
    "\n",
    "### Environment setup\n",
    "\n",
    "This notebook works with or without an API key, but you will be given more permissive access to the data with an API key.\n",
    "\n",
    "* The [Planetary Computer Hub](https://planetarycomputer.microsoft.com/compute) is pre-configured to use your API key.\n",
    "* To use your API key locally, set the environment variable `PC_SDK_SUBSCRIPTION_KEY` or use `planetary_computer.settings.set_subscription_key(<YOUR API Key>)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f3f030fa-8639-4b4e-b0a1-f3226375c28e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pystac_client\n",
    "from pystac.extensions.item_assets import ItemAssetsExtension\n",
    "import planetary_computer\n",
    "\n",
    "import matplotlib.colors\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import requests\n",
    "import rioxarray\n",
    "import rich.table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "527d242e-32b4-49fb-a513-a74ee38495da",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Query for data of interest\n",
    "\n",
    "We'll query for 2013 Chesapeake Land Use data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7e037147-5996-4532-a02c-5ce3d15b9a68",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fetching 2013\n",
      "Returned 166 Items\n"
     ]
    }
   ],
   "source": [
    "catalog = pystac_client.Client.open(\n",
    "    \"https://planetarycomputer.microsoft.com/api/stac/v1/\",\n",
    "    modifier=planetary_computer.sign_inplace,\n",
    ")\n",
    "\n",
    "latitude = 39.29\n",
    "longitude = -76.61\n",
    "datetimes = [\"2013\"]\n",
    "\n",
    "buffer = 0.6\n",
    "bbox = [longitude - buffer, latitude - buffer, longitude + buffer, latitude + buffer]\n",
    "items = dict()\n",
    "\n",
    "for datetime in datetimes:\n",
    "    print(f\"Fetching {datetime}\")\n",
    "    search = catalog.search(collections=[\"chesapeake-lu\"], bbox=bbox, datetime=datetime)\n",
    "\n",
    "items = list(search.get_items())\n",
    "print(f\"Returned {len(items)} Items\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49363704-3dd9-4f50-a8d2-b4b4c61f8360",
   "metadata": {},
   "source": [
    "Each Item contains a single \"data\" asset with a URL to the location of the Asset data on [Azure Blob Storage](https://azure.microsoft.com/en-us/services/storage/blobs/). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7bd73779-0545-4368-bf96-3d9cce489f97",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://landcoverdata.blob.core.windows.net/chesapeake/lu/BayWide_1m_LU_E1700000_N1960000.tif?st=2022-08-08T16%3A37%3A59Z&se=2022-08-16T16%3A37%3A59Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2022-08-09T16%3A37%3A58Z&ske=2022-08-16T16%3A37%3A58Z&sks=b&skv=2021-06-08&sig=8UW6710XVMKmP7BvKT4ZuULp87IR0bct2suXgk98RHA%3D\n"
     ]
    }
   ],
   "source": [
    "asset_href = items[0].assets[\"data\"].href\n",
    "print(asset_href)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6ccf07a-33bd-4be7-85c1-4471a17ae176",
   "metadata": {},
   "source": [
    "### Available STAC assets and metadata\n",
    "\n",
    "Let's check the available assets and metadata for the Chesapeake Land Use item. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a3f2152a-6925-4a01-b24f-485cca0b9f55",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃<span style=\"font-weight: bold\"> Asset key        </span>┃<span style=\"font-weight: bold\"> Title                           </span>┃\n",
       "┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ data             │ Land Use Classes                │\n",
       "│ tilejson         │ TileJSON with default rendering │\n",
       "│ rendered_preview │ Rendered preview                │\n",
       "└──────────────────┴─────────────────────────────────┘\n",
       "</pre>\n"
      ],
      "text/plain": [
       "┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃\u001b[1m \u001b[0m\u001b[1mAsset key       \u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mTitle                          \u001b[0m\u001b[1m \u001b[0m┃\n",
       "┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ data             │ Land Use Classes                │\n",
       "│ tilejson         │ TileJSON with default rendering │\n",
       "│ rendered_preview │ Rendered preview                │\n",
       "└──────────────────┴─────────────────────────────────┘\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = rich.table.Table(\"Asset key\", \"Title\")\n",
    "for k, asset in items[0].assets.items():\n",
    "    t.add_row(k, asset.title)\n",
    "t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "485ef558-d85e-41b4-898b-d799f205d45f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">┏━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃<span style=\"font-weight: bold\"> Key            </span>┃<span style=\"font-weight: bold\"> Value                                                                                          </span>┃\n",
       "┡━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ created        │ 2022-06-01T20:15:13.334677Z                                                                    │\n",
       "│ datetime       │ None                                                                                           │\n",
       "│ end_datetime   │ 2014-12-31T23:59:59Z                                                                           │\n",
       "│ proj:epsg      │ None                                                                                           │\n",
       "│ proj:shape     │ [10000, 10000]                                                                                 │\n",
       "│ proj:transform │ [1.0, 0.0, 1699999.1624823206, 0.0, -1.0, 1970000.8902923968]                                  │\n",
       "│ proj:wkt2      │ PROJCS[\"USA_Contiguous_Albers_Equal_Area_Conic_USGS_version\",GEOGCS[\"NAD83\",DATUM[\"North_Amer… │\n",
       "│                │ 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Gr… │\n",
       "│ start_datetime │ 2013-01-01T00:00:00Z                                                                           │\n",
       "└────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────┘\n",
       "</pre>\n"
      ],
      "text/plain": [
       "┏━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃\u001b[1m \u001b[0m\u001b[1mKey           \u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mValue                                                                                         \u001b[0m\u001b[1m \u001b[0m┃\n",
       "┡━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ created        │ 2022-06-01T20:15:13.334677Z                                                                    │\n",
       "│ datetime       │ None                                                                                           │\n",
       "│ end_datetime   │ 2014-12-31T23:59:59Z                                                                           │\n",
       "│ proj:epsg      │ None                                                                                           │\n",
       "│ proj:shape     │ [10000, 10000]                                                                                 │\n",
       "│ proj:transform │ [1.0, 0.0, 1699999.1624823206, 0.0, -1.0, 1970000.8902923968]                                  │\n",
       "│ proj:wkt2      │ PROJCS[\"USA_Contiguous_Albers_Equal_Area_Conic_USGS_version\",GEOGCS[\"NAD83\",DATUM[\"North_Amer… │\n",
       "│                │ 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Gr… │\n",
       "│ start_datetime │ 2013-01-01T00:00:00Z                                                                           │\n",
       "└────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────┘\n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = rich.table.Table(\"Key\", \"Value\")\n",
    "for k, v in sorted(items[0].properties.items()):\n",
    "    t.add_row(k, str(v))\n",
    "\n",
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f76a3b43-7827-4e8d-9305-6af33794bbce",
   "metadata": {},
   "source": [
    "### Displaying the data\n",
    "\n",
    "This dataset includes a preferred colormap mapping raster values to colors. The Collection's `item_assets` field includes an overview of the class descriptions and values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cd296765-a455-4a92-bd1f-4e29da0344ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━┓\n",
       "┃<span style=\"font-weight: bold\"> Description                          </span>┃<span style=\"font-weight: bold\"> Value </span>┃\n",
       "┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━┩\n",
       "│ Impervious Roads                     │ 1     │\n",
       "│ Impervious Non-Roads                 │ 2     │\n",
       "│ Tree Canopy over Impervious Surfaces │ 3     │\n",
       "│ Water                                │ 4     │\n",
       "│ Tidal Wetlands                       │ 5     │\n",
       "│ Floodplain Wetlands                  │ 6     │\n",
       "│ Other Wetlands                       │ 7     │\n",
       "│ Forest                               │ 8     │\n",
       "│ Tree Canopy over Turf Grass          │ 9     │\n",
       "│ Mixed Open                           │ 10    │\n",
       "│ Fractional Turf (small)              │ 11    │\n",
       "│ Fractional Turf (medium)             │ 12    │\n",
       "│ Fractional Turf (large)              │ 13    │\n",
       "│ Fractional Impervious                │ 14    │\n",
       "│ Turf Grass                           │ 15    │\n",
       "│ Cropland                             │ 16    │\n",
       "│ Pasture/Hay                          │ 17    │\n",
       "└──────────────────────────────────────┴───────┘\n",
       "</pre>\n"
      ],
      "text/plain": [
       "┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━┓\n",
       "┃\u001b[1m \u001b[0m\u001b[1mDescription                         \u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mValue\u001b[0m\u001b[1m \u001b[0m┃\n",
       "┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━┩\n",
       "│ Impervious Roads                     │ 1     │\n",
       "│ Impervious Non-Roads                 │ 2     │\n",
       "│ Tree Canopy over Impervious Surfaces │ 3     │\n",
       "│ Water                                │ 4     │\n",
       "│ Tidal Wetlands                       │ 5     │\n",
       "│ Floodplain Wetlands                  │ 6     │\n",
       "│ Other Wetlands                       │ 7     │\n",
       "│ Forest                               │ 8     │\n",
       "│ Tree Canopy over Turf Grass          │ 9     │\n",
       "│ Mixed Open                           │ 10    │\n",
       "│ Fractional Turf (small)              │ 11    │\n",
       "│ Fractional Turf (medium)             │ 12    │\n",
       "│ Fractional Turf (large)              │ 13    │\n",
       "│ Fractional Impervious                │ 14    │\n",
       "│ Turf Grass                           │ 15    │\n",
       "│ Cropland                             │ 16    │\n",
       "│ Pasture/Hay                          │ 17    │\n",
       "└──────────────────────────────────────┴───────┘\n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Classifications\n",
    "collection = catalog.get_collection(\"chesapeake-lu\")\n",
    "ia = ItemAssetsExtension.ext(collection)\n",
    "x = ia.item_assets[\"data\"]\n",
    "\n",
    "class_names = {\n",
    "    x[\"description\"]: x[\"value\"] for x in x.properties[\"classification:classes\"]\n",
    "}\n",
    "class_values = {v: k for k, v in class_names.items()}\n",
    "\n",
    "t = rich.table.Table(\"Description\", \"Value\")\n",
    "for k, v in class_names.items():\n",
    "    t.add_row(k, str(v))\n",
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e8867d7-e6a2-447a-970f-b1a682239655",
   "metadata": {},
   "source": [
    "And the Planetary Computer's [Data API](https://planetarycomputer.microsoft.com/api/data/v1/docs) includes the colormap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "55f58de8-36e4-4db8-98e8-d69f3a2de749",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'0': [0, 0, 0, 0],\n",
       " '1': [0, 0, 0, 255],\n",
       " '2': [115, 0, 0, 255],\n",
       " '3': [85, 255, 0, 255],\n",
       " '4': [0, 112, 255, 255],\n",
       " '5': [0, 255, 197, 255],\n",
       " '6': [0, 230, 169, 255],\n",
       " '7': [0, 230, 169, 255],\n",
       " '8': [38, 115, 0, 255],\n",
       " '9': [170, 255, 0, 255],\n",
       " '10': [168, 112, 0, 255],\n",
       " '11': [255, 190, 232, 255],\n",
       " '12': [255, 190, 232, 255],\n",
       " '13': [255, 190, 232, 255],\n",
       " '14': [197, 0, 255, 255],\n",
       " '15': [255, 255, 115, 255],\n",
       " '16': [230, 152, 0, 255],\n",
       " '17': [230, 152, 0, 255],\n",
       " '255': [0, 0, 0, 0]}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "classmap = requests.get(\n",
    "    \"https://pct-apis-staging.westeurope.cloudapp.azure.com/data/legend/classmap/chesapeake-lu\"  # noqa: E501\n",
    ").json()\n",
    "classmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cae2e2d4-9c5b-46e7-8f18-7bd78fee055b",
   "metadata": {},
   "source": [
    "We'll convert those values to a [matplotlib Colormap](https://matplotlib.org/stable/api/_as_gen/matplotlib.colors.ListedColormap.html) for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ef1d08dd-0427-4489-a512-7a0d99da5e97",
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = [matplotlib.colors.to_rgba([x / 255 for x in c]) for c in classmap.values()]\n",
    "cmap = matplotlib.colors.ListedColormap(colors, name=\"chesapeak-lu\")\n",
    "ticks = np.arange(cmap.N)\n",
    "labels = [class_values.get(value, \"nodata\") for value in ticks]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4a4cab6-b221-4f5f-9536-dc80d6ce86b2",
   "metadata": {},
   "source": [
    "Finally, we can read the data and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "82426bf9-db88-45f1-963a-e671e6751d56",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = rioxarray.open_rasterio(items[0].assets[\"data\"].href).squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "47c95711-fc9f-4591-bfdd-1b93b6521c89",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-chesapeake-chesapeake-example.ipynb/10.png\"/>"
      ],
      "text/plain": [
       "<Figure size 1152x864 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 12))\n",
    "\n",
    "ds.isel(x=slice(3000), y=slice(3000)).plot(\n",
    "    ax=ax, cmap=cmap, norm=matplotlib.colors.NoNorm()\n",
    ")\n",
    "ax.set_axis_off()\n",
    "ax.set(title=None)\n",
    "\n",
    "colorbar = fig.axes[1]\n",
    "colorbar.set_yticks(ticks, labels=labels);"
   ]
  }
 ],
 "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.8.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}