{ "cells": [ { "cell_type": "markdown", "id": "d2583da9-3f1e-40e6-b81d-700f26e89daa", "metadata": { "tags": [] }, "source": [ "## Accessing NOAA C-CAP Regional Land Cover with the Planetary Computer STAC API\n", "\n", "The [National Oceanic and Atmospheric Administration (NOAA) Coastal Change Analysis Program (C-CAP) Regional Land Cover](https://coast.noaa.gov/digitalcoast/data/ccapregional.html) dataset features standardized raster based land cover for coastal areas of the contiguous United States (CONUS) and Hawaii at a 30-meter resolution available from 1975 through 2016. The data are derived from Landsat's Thematic Mapper (TM) satellite imagery and ancillary information.\n", "\n", "The use of standardized data and procedures assures consistency through time and across geographies. C-CAP data forms the coastal expression of the National Land Cover Database (NLCD) and the A-16 land cover theme of the National Spatial Data Infrastructure (NSDI). \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/noaa-c-cap)." ] }, { "cell_type": "markdown", "id": "6c7da081-b792-45f2-a428-ea8f9ecbda08", "metadata": {}, "source": [ "### 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", "The [Planetary Computer Hub](https://planetarycomputer.microsoft.com/compute) is pre-configured to use your API key." ] }, { "cell_type": "code", "execution_count": 1, "id": "f3f030fa-8639-4b4e-b0a1-f3226375c28e", "metadata": {}, "outputs": [], "source": [ "from matplotlib.colors import ListedColormap\n", "from pystac.extensions.item_assets import ItemAssetsExtension\n", "\n", "import numpy as np\n", "import odc.stac\n", "import planetary_computer\n", "import pystac_client\n", "import rasterio\n", "import rasterio.features\n", "import rich.table" ] }, { "cell_type": "markdown", "id": "980b87ca-9fcb-4146-94ef-420542debf7c", "metadata": {}, "source": [ "### Data access\n", "\n", "The datasets hosted by the Planetary Computer are available from [Azure Blob Storage](https://docs.microsoft.com/en-us/azure/storage/blobs/). We'll use [pystac-client](https://pystac-client.readthedocs.io/) to search the Planetary Computer's [STAC API](https://planetarycomputer.microsoft.com/api/stac/v1/docs) for the subset of the data that we care about, and then we'll load the data directly from Azure Blob Storage. We'll specify a `modifier` so that we can access the data stored in the Planetary Computer's private Blob Storage Containers. See [Reading from the STAC API](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/) and [Using tokens for data access](https://planetarycomputer.microsoft.com/docs/concepts/sas/) for more." ] }, { "cell_type": "code", "execution_count": 2, "id": "66b5faa1-abbe-4c3f-8b86-90b83722b9b0", "metadata": {}, "outputs": [], "source": [ "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")" ] }, { "cell_type": "markdown", "id": "527d242e-32b4-49fb-a513-a74ee38495da", "metadata": { "tags": [] }, "source": [ "### Query for available data\n", "\n", "The NOAA C-CAP dataset reaches back to 1975 in select areas. We selected the coordinates for Lansing, MI to showcase data for 1975 and use the STAC API to find what data items are available." ] }, { "cell_type": "code", "execution_count": 3, "id": "63ddd01e-2ae3-4b5d-af8f-cdcbb56f0423", "metadata": {}, "outputs": [], "source": [ "# Select area and time of interest\n", "latitude = 42.7325\n", "longitude = -84.5555\n", "location = [longitude, latitude]\n", "geometry = {\n", " \"type\": \"Point\",\n", " \"coordinates\": location,\n", "}\n", "datetimes = [\n", " \"1975\",\n", "]\n", "\n", "buffer = 2\n", "bbox = [longitude - buffer, latitude - buffer, longitude + buffer, latitude + buffer]\n", "items = dict()\n", "search = catalog.search(collections=\"noaa-c-cap\", intersects=geometry, datetime=\"1975\")\n", "item = next(search.items())" ] }, { "cell_type": "markdown", "id": "69f3564c-bd52-452d-875d-34d9309a64c9", "metadata": {}, "source": [ "### Available Assets" ] }, { "cell_type": "markdown", "id": "cfd0f626-55bb-4bcf-a2b6-d276df55a0d4", "metadata": {}, "source": [ "Let's display the available assets for the NOAA C-CAP item. " ] }, { "cell_type": "code", "execution_count": 4, "id": "a3f2152a-6925-4a01-b24f-485cca0b9f55", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃ Key              ┃ Title                           ┃\n",
       "┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ data             │                                 │\n",
       "│ tilejson         │ TileJSON with default rendering │\n",
       "│ rendered_preview │ Rendered preview                │\n",
       "└──────────────────┴─────────────────────────────────┘\n",
       "
\n" ], "text/plain": [ "┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n", "┃\u001b[1m \u001b[0m\u001b[1mKey \u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mTitle \u001b[0m\u001b[1m \u001b[0m┃\n", "┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n", "│ data │ │\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(\"Key\", \"Title\")\n", "for key, asset in item.assets.items():\n", " t.add_row(key, asset.title)\n", "t" ] }, { "cell_type": "markdown", "id": "fcd49bc9-32ac-4b10-83f8-e9fbd2a870ac", "metadata": {}, "source": [ "### Loading the land cover data\n", "\n", "For this example, we'll visualize the land cover classifications surrounding Lansing, MI. Let's grab the data COG and load them into an xarray using [odc-stac](https://github.com/opendatacube/odc-stac). \n" ] }, { "cell_type": "code", "execution_count": 5, "id": "a8fc2966-2787-4f20-bb2d-e47d4829cec5", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:      (y: 16091, x: 12565, time: 1)\n",
       "Coordinates:\n",
       "  * y            (y) float64 2.491e+06 2.491e+06 ... 2.008e+06 2.008e+06\n",
       "  * x            (x) float64 7.454e+05 7.455e+05 ... 1.122e+06 1.122e+06\n",
       "    spatial_ref  int32 5070\n",
       "  * time         (time) datetime64[ns] 1975-01-01\n",
       "Data variables:\n",
       "    data         (time, y, x) uint8 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
" ], "text/plain": [ "\n", "Dimensions: (y: 16091, x: 12565, time: 1)\n", "Coordinates:\n", " * y (y) float64 2.491e+06 2.491e+06 ... 2.008e+06 2.008e+06\n", " * x (x) float64 7.454e+05 7.455e+05 ... 1.122e+06 1.122e+06\n", " spatial_ref int32 5070\n", " * time (time) datetime64[ns] 1975-01-01\n", "Data variables:\n", " data (time, y, x) uint8 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = odc.stac.load(\n", " [item],\n", " crs=\"EPSG:5070\",\n", " bbox=bbox,\n", " bands=\"data\",\n", " resolution=30,\n", ")\n", "data" ] }, { "cell_type": "markdown", "id": "f76a3b43-7827-4e8d-9305-6af33794bbce", "metadata": {}, "source": [ "### Displaying the data\n", "\n", "Let's use the classification colors as defined in the [C-CAP Classification Scheme and Class Definitions](https://coast.noaa.gov/data/digitalcoast/pdf/ccap-class-scheme-regional.pdf) document. The classification colors are also defined in a colormap contained in the source GeoTIFF image. We'll extract the colormap from the GeoTIFF image and apply it to the plot." ] }, { "cell_type": "code", "execution_count": 6, "id": "cd296765-a455-4a92-bd1f-4e29da0344ab", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'Background': 0,\n", " 'Unclassified (Cloud, Shadow, etc)': 1,\n", " 'High Intensity Developed': 2,\n", " 'Medium Intensity Developed': 3,\n", " 'Low Intensity Developed': 4,\n", " 'Developed Open Space': 5,\n", " 'Cultivated Crops': 6,\n", " 'Pasture/Hay': 7,\n", " 'Grassland/Herbaceous': 8,\n", " 'Deciduous Forest': 9,\n", " 'Evergreen Forest': 10,\n", " 'Mixed Forest': 11,\n", " 'Scrub/Shrub': 12,\n", " 'Palustrine Forested Wetland': 13,\n", " 'Palustrine Scrub/Shrub Wetland': 14,\n", " 'Palustrine Emergent Wetland (Persistent)': 15,\n", " 'Estuarine Forested Wetland': 16,\n", " 'Estuarine Scrub/Shrub Wetland': 17,\n", " 'Estuarine Emergent Wetland': 18,\n", " 'Unconsolidated Shore': 19,\n", " 'Bare Land': 20,\n", " 'Open Water': 21,\n", " 'Palustrine Aquatic Bed': 22,\n", " 'Estuarine Aquatic Bed': 23,\n", " 'Tundra': 24,\n", " 'Perennial Ice/Snow': 25}" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Extract colormap from GeoTIFF\n", "collection = catalog.get_collection(\"noaa-c-cap\")\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", "values_to_classes = {v: k for k, v in class_names.items()}\n", "class_count = len(class_names)\n", "class_names" ] }, { "cell_type": "code", "execution_count": 7, "id": "533ff2aa-bd4a-4fb5-a0bb-1b1118a7cb8a", "metadata": {}, "outputs": [], "source": [ "# Define colormap from STAC data\n", "with rasterio.open(item.assets[\"data\"].href) as src:\n", " colormap_def = src.colormap(1) # get metadata colormap for band 1\n", " colormap = [\n", " np.array(colormap_def[i]) / 255 for i in range(len(colormap_def))\n", " ] # transform to matplotlib color format\n", "\n", "cmap = ListedColormap(colormap)" ] }, { "cell_type": "code", "execution_count": 8, "id": "e83c096f-0ef5-46a2-98e5-2b59d116004b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Visualization\n", "g = data[\"data\"].plot.imshow(\n", " cmap=cmap,\n", " col=\"time\",\n", " vmin=0,\n", " vmax=(len(colormap_def) - 1),\n", " col_wrap=1,\n", " size=8,\n", ")\n", "datetimes = data[\"data\"].time.to_pandas().dt.strftime(\"%Y\")\n", "\n", "for ax, datetime in zip(g.axes.flat, datetimes):\n", " ax.set_title(datetime)" ] }, { "cell_type": "markdown", "id": "f10b7df2-0064-43d7-b8d6-1862caaa6b36", "metadata": {}, "source": [ "Rendering notes:\n", "1. The `vmin` and `vmax` values are set to match the number of classes defined in the colormap in the source GeoTIFF image. Even though there are only 26 classes defined in the [C-CAP Classification Scheme and Class Definitions](https://coast.noaa.gov/data/digitalcoast/pdf/ccap-class-scheme-regional.pdf) document, the GeoTIFF produces 256 class values and corresponding colors. " ] } ], "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 }