{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Accessing MODIS snow data with the Planetary Computer STAC API\n", "\n", "The planetary computer hosts two snow-related MODIS 6.1 products:\n", "\n", "- Snow cover daily (10A1)\n", "- Snow cover 8-day (10A2)\n", "\n", "For more information about the products themselves, check out the Data Pages at the [bottom of this document](#data-pages)." ] }, { "cell_type": "markdown", "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 is pre-configured to use your API key." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import odc.stac\n", "import planetary_computer\n", "import pystac_client\n", "import rich.table" ] }, { "cell_type": "markdown", "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, "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", "metadata": {}, "source": [ "### Query for available data\n", "\n", "MODIS is a global dataset with a variety of products available within each larger category (vegetation, snow, fire, temperature, and reflectance).\n", "The [MODIS group](https://planetarycomputer.microsoft.com/dataset/group/modis) contains a complete listing of available collections.\n", "Each collection's id is in the format `modis-{product}-061`, where `product` is the MODIS product id.\n", "The `-061` suffix indicates that all of the MODIS collections are part of the [MODIS 6.1 update](https://atmosphere-imager.gsfc.nasa.gov/documentation/collection-61).\n", "\n", "We'll look at the snow cover around the Mount of the Holy Cross in Colorado and how it progresses throughout the winter, using the daily snow product (10A1)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fetching 2020-11\n", "Fetching 2020-12\n", "Fetching 2021-01\n", "Fetching 2021-02\n", "Fetching 2021-03\n", "Fetching 2021-04\n", "Fetching 2021-05\n", "Fetching 2021-06\n", "{'2020-11': , '2020-12': , '2021-01': , '2021-02': , '2021-03': , '2021-04': , '2021-05': , '2021-06': }\n" ] } ], "source": [ "longitude = -106.481687\n", "latitude = 39.466829\n", "mount_of_the_holy_cross = [longitude, latitude]\n", "geometry = {\n", " \"type\": \"Point\",\n", " \"coordinates\": mount_of_the_holy_cross,\n", "}\n", "datetimes = [\n", " \"2020-11\",\n", " \"2020-12\",\n", " \"2021-01\",\n", " \"2021-02\",\n", " \"2021-03\",\n", " \"2021-04\",\n", " \"2021-05\",\n", " \"2021-06\",\n", "]\n", "items = dict()\n", "\n", "for datetime in datetimes:\n", " print(f\"Fetching {datetime}\")\n", " search = catalog.search(\n", " collections=[\"modis-10A1-061\"],\n", " intersects=geometry,\n", " datetime=datetime,\n", " )\n", " items[datetime] = search.get_all_items()[0]\n", "\n", "print(items)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Available assets\n", "\n", "Each item has several available assets, including the original HDF file and a Cloud-optimized GeoTIFF of each subdataset." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃ Key                                ┃ Title                                                                   ┃\n",
       "┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ hdf                                │ Source data containing all bands                                        │\n",
       "│ NDSI                               │ Raw NDSI values (i.e. prior to screening).                              │\n",
       "│ metadata                           │ Federal Geographic Data Committee (FGDC) Metadata                       │\n",
       "│ orbit_pnt                          │ Pointer to the orbit of the swath mapped into each grid cell.           │\n",
       "│ granule_pnt                        │ Pointer for identifying the swath mapped into each grid cell.           │\n",
       "│ NDSI_Snow_Cover                    │ Gridded NDSI snow cover and data flag values.                           │\n",
       "│ Snow_Albedo_Daily_Tile             │ Daily snow albedo corresponding to the NDSI_Snow_Cover parameter.       │\n",
       "│ NDSI_Snow_Cover_Basic_QA           │ A general estimate of the quality of the algorithm result.              │\n",
       "│ NDSI_Snow_Cover_Algorithm_Flags_QA │ Algorithm-specific bit flags set for data screens and for inland water. │\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", "│ hdf │ Source data containing all bands │\n", "│ NDSI │ Raw NDSI values (i.e. prior to screening). │\n", "│ metadata │ Federal Geographic Data Committee (FGDC) Metadata │\n", "│ orbit_pnt │ Pointer to the orbit of the swath mapped into each grid cell. │\n", "│ granule_pnt │ Pointer for identifying the swath mapped into each grid cell. │\n", "│ NDSI_Snow_Cover │ Gridded NDSI snow cover and data flag values. │\n", "│ Snow_Albedo_Daily_Tile │ Daily snow albedo corresponding to the NDSI_Snow_Cover parameter. │\n", "│ NDSI_Snow_Cover_Basic_QA │ A general estimate of the quality of the algorithm result. │\n", "│ NDSI_Snow_Cover_Algorithm_Flags_QA │ Algorithm-specific bit flags set for data screens and for inland water. │\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 items[\"2020-11\"].assets.items():\n", " t.add_row(key, asset.title)\n", "t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Loading the snow cover data\n", "\n", "For this example, we'll visualize and compare the snow cover month to month.\n", "Let's grab each snow cover COG and load them into an xarray using [odc-stac](https://github.com/opendatacube/odc-stac).\n", "We'll also apply the scaling as defined by the `raster:bands` extension.\n", "The MODIS coordinate reference system is a [sinusoidal grid](https://modis-land.gsfc.nasa.gov/MODLAND_grid.html), which means that views in a naïve XY raster look skewed.\n", "For visualization purposes, we reproject to a [spherical Mercator projection](https://wiki.openstreetmap.org/wiki/EPSG:3857) for intuitive, north-up visualization.\n", "\n", "The NDSI Snow Cover values are defined as:\n", "\n", "```\n", "0–100: NDSI snow cover\n", "200: missing data\n", "201: no decision\n", "211: night\n", "237: inland water\n", "239: ocean\n", "250: cloud\n", "254: detector saturated\n", "255: fill\n", "```\n", "\n", "We want to mask out all numbers greater than 100." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:          (time: 8, y: 290, x: 224)\n",
       "Coordinates:\n",
       "  * y                (y) float64 4.861e+06 4.861e+06 ... 4.717e+06 4.717e+06\n",
       "  * x                (x) float64 -1.191e+07 -1.191e+07 ... -1.18e+07 -1.18e+07\n",
       "    spatial_ref      int32 3857\n",
       "  * time             (time) datetime64[ns] 2020-11-30 2020-12-31 ... 2021-06-30\n",
       "Data variables:\n",
       "    NDSI_Snow_Cover  (time, y, x) float64 nan nan nan nan ... nan nan nan nan
" ], "text/plain": [ "\n", "Dimensions: (time: 8, y: 290, x: 224)\n", "Coordinates:\n", " * y (y) float64 4.861e+06 4.861e+06 ... 4.717e+06 4.717e+06\n", " * x (x) float64 -1.191e+07 -1.191e+07 ... -1.18e+07 -1.18e+07\n", " spatial_ref int32 3857\n", " * time (time) datetime64[ns] 2020-11-30 2020-12-31 ... 2021-06-30\n", "Data variables:\n", " NDSI_Snow_Cover (time, y, x) float64 nan nan nan nan ... nan nan nan nan" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bbox = [longitude - 0.5, latitude - 0.5, longitude + 0.5, latitude + 0.5]\n", "data = odc.stac.load(\n", " items.values(),\n", " crs=\"EPSG:3857\",\n", " bbox=bbox,\n", " bands=\"NDSI_Snow_Cover\",\n", " resolution=500,\n", ")\n", "data = data.where(data <= 100, drop=True)\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Displaying the data\n", "\n", "Let's display the snow cover for each month." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "g = data[\"NDSI_Snow_Cover\"].plot.imshow(\n", " col=\"time\", vmin=0, vmax=100, col_wrap=4, size=4\n", ")\n", "months = data[\"NDSI_Snow_Cover\"].time.to_pandas().dt.strftime(\"%B %Y\")\n", "\n", "for ax, month in zip(g.axes.flat, months):\n", " ax.set_title(month)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You'll notice there's a lot of missing values due to masking, probably due to clouds.\n", "More sophisticated analysis would use the `NDSI_Snow_Cover_Basic_QA` and other bands to merge information from multiple scenes, or pick the best scenes for analysis.\n", "\n", "### Data pages\n", "\n", "These pages include links to the user guides:\n", "\n", "- MOD10A1: https://nsidc.org/data/MOD10A1/versions/61\n", "- MOD10A2: https://nsidc.org/data/MOD10A2/versions/61" ] } ], "metadata": { "interpreter": { "hash": "33388988d2b1cb540819d28e8f67ee9e7e7397cf80503e5dd9c105f19f832c00" }, "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": 4 }