{ "cells": [ { "cell_type": "markdown", "id": "0d97ffaf-9b10-423c-a505-e0b19949844e", "metadata": {}, "source": [ "## Accessing Landsat Collection 2 Level-1 and Level-2 data with the Planetary Computer STAC API\n", "\n", "The Planetary Computer's [Landsat](https://landsat.gsfc.nasa.gov/) dataset represents a global archive of [Level-1](https://www.usgs.gov/landsat-missions/landsat-collection-2-level-1-data) and [Level-2](https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products) data from from [Landsat Collection 2](https://www.usgs.gov/core-science-systems/nli/landsat/landsat-collection-2). The dataset is grouped into two STAC Collections:\n", "- `landsat-c2-l2` for Level-2 data\n", "- `landsat-c2-l1` for Level-1 data\n", "\n", "This notebook demonstrates the use of the Planetary Computer STAC API to query for Landsat data. We will start with an example using Level-2 data and follow with a similar example using Level-1 data.\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", "- 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 `pc.settings.set_subscription_key()`." ] }, { "cell_type": "code", "execution_count": 1, "id": "f50f9f62", "metadata": {}, "outputs": [], "source": [ "import pystac_client\n", "import planetary_computer\n", "import odc.stac\n", "import matplotlib.pyplot as plt\n", "\n", "from pystac.extensions.eo import EOExtension as eo" ] }, { "cell_type": "markdown", "id": "4e154dc2-ae73-4ca5-9918-7153aea89919", "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": "bfe8e3d9-035c-4928-8752-e0329d87a7b0", "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": "198b54e7", "metadata": {}, "source": [ "### Choose an area and time of interest\n", "\n", "This area is in Redmond, WA, USA, near Microsoft's main campus. We'll search all of 2021, limiting the results to those less than 10% cloudy." ] }, { "cell_type": "code", "execution_count": 3, "id": "d1a949e6", "metadata": {}, "outputs": [], "source": [ "bbox_of_interest = [-122.2751, 47.5469, -121.9613, 47.7458]\n", "time_of_interest = \"2021-01-01/2021-12-31\"" ] }, { "cell_type": "code", "execution_count": 4, "id": "a56a65b0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Returned 13 Items\n" ] } ], "source": [ "search = catalog.search(\n", " collections=[\"landsat-c2-l2\"],\n", " bbox=bbox_of_interest,\n", " datetime=time_of_interest,\n", " query={\"eo:cloud_cover\": {\"lt\": 10}},\n", ")\n", "\n", "items = search.item_collection()\n", "print(f\"Returned {len(items)} Items\")" ] }, { "cell_type": "markdown", "id": "f0a8210f", "metadata": {}, "source": [ "Let's find the least cloudy of the bunch." ] }, { "cell_type": "code", "execution_count": 5, "id": "6204ef0c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Choosing LC08_L2SP_046027_20210725_02_T1 from 2021-07-25 with 0.38% cloud cover\n" ] } ], "source": [ "selected_item = min(items, key=lambda item: eo.ext(item).cloud_cover)\n", "\n", "print(\n", " f\"Choosing {selected_item.id} from {selected_item.datetime.date()}\"\n", " + f\" with {selected_item.properties['eo:cloud_cover']}% cloud cover\"\n", ")" ] }, { "cell_type": "markdown", "id": "e2a89114", "metadata": {}, "source": [ "### Available assets\n", "\n", "In addition to numerous metadata assets, each Electro-Optical (EO) band is a separate asset." ] }, { "cell_type": "code", "execution_count": 6, "id": "272bdccc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " qa: Surface Temperature Quality Assessment Band\n", " ang: Angle Coefficients File\n", " red: Red Band\n", " blue: Blue Band\n", " drad: Downwelled Radiance Band\n", " emis: Emissivity Band\n", " emsd: Emissivity Standard Deviation Band\n", " trad: Thermal Radiance Band\n", " urad: Upwelled Radiance Band\n", " atran: Atmospheric Transmittance Band\n", " cdist: Cloud Distance Band\n", " green: Green Band\n", " nir08: Near Infrared Band 0.8\n", " lwir11: Surface Temperature Band\n", " swir16: Short-wave Infrared Band 1.6\n", " swir22: Short-wave Infrared Band 2.2\n", " coastal: Coastal/Aerosol Band\n", " mtl.txt: Product Metadata File (txt)\n", " mtl.xml: Product Metadata File (xml)\n", " mtl.json: Product Metadata File (json)\n", " qa_pixel: Pixel Quality Assessment Band\n", " qa_radsat: Radiometric Saturation and Terrain Occlusion Quality Assessment Band\n", " qa_aerosol: Aerosol Quality Assessment Band\n", " tilejson: TileJSON with default rendering\n", "rendered_preview: Rendered preview\n" ] } ], "source": [ "max_key_length = len(max(selected_item.assets, key=len))\n", "for key, asset in selected_item.assets.items():\n", " print(f\"{key.rjust(max_key_length)}: {asset.title}\")" ] }, { "cell_type": "markdown", "id": "2001fcf4", "metadata": {}, "source": [ "### Render a natural color image of the AOI\n", "\n", "We'll start by loading the red, green, and blue bands for our area of interest into an xarray dataset using [`odc-stac`](https://pypi.org/project/odc-stac/). We will also load the nir08 band for use in computing an NDVI value in a later example. Note that we pass the `sign` function from the [planetary-computer](https://github.com/microsoft/planetary-computer-sdk-for-python) package so that `odc-stac` can supply the required [Shared Access Signature](https://docs.microsoft.com/en-us/azure/storage/common/storage-sas-overview) tokens that are necessary to download the asset data." ] }, { "cell_type": "code", "execution_count": 7, "id": "9849f2e0", "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: 747, x: 795)\n",
       "Coordinates:\n",
       "  * y            (y) float64 5.289e+06 5.289e+06 ... 5.266e+06 5.266e+06\n",
       "  * x            (x) float64 5.543e+05 5.544e+05 ... 5.781e+05 5.781e+05\n",
       "    spatial_ref  int32 32610\n",
       "    time         datetime64[ns] 2021-07-25T18:55:39.475647\n",
       "Data variables:\n",
       "    nir08        (y, x) uint16 7191 7187 7216 7237 ... 19509 21532 20871 20917\n",
       "    red          (y, x) uint16 7172 7158 7193 7219 7246 ... 8067 8171 9372 8634\n",
       "    green        (y, x) uint16 7498 7489 7518 7541 7570 ... 8706 8771 9580 9062\n",
       "    blue         (y, x) uint16 7337 7343 7370 7433 7467 ... 8011 8068 8616 8252\n",
       "    qa_pixel     (y, x) uint16 21952 21952 21952 21952 ... 21824 21824 21824\n",
       "    lwir11       (y, x) uint16 43477 43466 43464 43468 ... 44336 44431 44355
" ], "text/plain": [ "\n", "Dimensions: (y: 747, x: 795)\n", "Coordinates:\n", " * y (y) float64 5.289e+06 5.289e+06 ... 5.266e+06 5.266e+06\n", " * x (x) float64 5.543e+05 5.544e+05 ... 5.781e+05 5.781e+05\n", " spatial_ref int32 32610\n", " time datetime64[ns] 2021-07-25T18:55:39.475647\n", "Data variables:\n", " nir08 (y, x) uint16 7191 7187 7216 7237 ... 19509 21532 20871 20917\n", " red (y, x) uint16 7172 7158 7193 7219 7246 ... 8067 8171 9372 8634\n", " green (y, x) uint16 7498 7489 7518 7541 7570 ... 8706 8771 9580 9062\n", " blue (y, x) uint16 7337 7343 7370 7433 7467 ... 8011 8068 8616 8252\n", " qa_pixel (y, x) uint16 21952 21952 21952 21952 ... 21824 21824 21824\n", " lwir11 (y, x) uint16 43477 43466 43464 43468 ... 44336 44431 44355" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bands_of_interest = [\"nir08\", \"red\", \"green\", \"blue\", \"qa_pixel\", \"lwir11\"]\n", "data = odc.stac.stac_load(\n", " [selected_item], bands=bands_of_interest, bbox=bbox_of_interest\n", ").isel(time=0)\n", "data" ] }, { "cell_type": "markdown", "id": "e0b24f78", "metadata": {}, "source": [ "Now we'll convert our xarray Dataset to a DataArray and plot the RGB image. We set `robust=True` in `imshow` to avoid manual computation of the color limits (vmin and vmax) that is necessary for data not scaled to between 0-1 while also eliminating extreme values that can cause a washed out image." ] }, { "cell_type": "code", "execution_count": 8, "id": "14686fde", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(10, 10))\n", "\n", "data[[\"red\", \"green\", \"blue\"]].to_array().plot.imshow(robust=True, ax=ax)\n", "ax.set_title(\"Natural Color, Redmond, WA\");" ] }, { "cell_type": "markdown", "id": "a9275433", "metadata": {}, "source": [ "### Render an NDVI image of the AOI\n", "\n", "Landsat has several bands, and with them we can go beyond rendering natural color imagery; for example, the following code computes a [Normalized Difference Vegetation Index (NDVI)](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index) using the near-infrared and red bands. Note that we convert the red and near infrared bands to a data type that can contain negative values; if this is not done, negative NDVI values will be incorrectly stored." ] }, { "cell_type": "code", "execution_count": 9, "id": "ffd2eefb", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "red = data[\"red\"].astype(\"float\")\n", "nir = data[\"nir08\"].astype(\"float\")\n", "ndvi = (nir - red) / (nir + red)\n", "\n", "fig, ax = plt.subplots(figsize=(14, 10))\n", "ndvi.plot.imshow(ax=ax, cmap=\"viridis\")\n", "ax.set_title(\"NDVI, Redmond, WA\");" ] }, { "cell_type": "markdown", "id": "e3d4704e-b825-418c-8809-957060623fc4", "metadata": {}, "source": [ "### Selecting specific platforms\n", "\n", "Landsat Collection 2 Level-2 includes assets from several different Platforms." ] }, { "cell_type": "code", "execution_count": 10, "id": "99da3fb2-12ec-47bc-94d8-99eb24298de5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['landsat-4', 'landsat-5', 'landsat-7', 'landsat-8', 'landsat-9']" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "catalog.get_collection(\"landsat-c2-l2\").summaries.to_dict()[\"platform\"]" ] }, { "cell_type": "markdown", "id": "686b3cf7-88fe-4baa-8513-8285ebf39eec", "metadata": {}, "source": [ "You might want to limit your search to a specific platform or platforms, to avoid the [Landsat 7 Scan Line Corrector failure](https://www.usgs.gov/landsat-missions/landsat-7), for example. Use the `\"platform\": {\"in\": ... }` query to select specific platforms." ] }, { "cell_type": "code", "execution_count": 11, "id": "b72663b7-2a5c-44fc-bca2-2c9dff6432f4", "metadata": {}, "outputs": [], "source": [ "search = catalog.search(\n", " collections=[\"landsat-c2-l2\"],\n", " bbox=bbox_of_interest,\n", " datetime=time_of_interest,\n", " query={\n", " \"eo:cloud_cover\": {\"lt\": 10},\n", " \"platform\": {\"in\": [\"landsat-8\", \"landsat-9\"]},\n", " },\n", ")\n", "items = search.get_all_items()" ] }, { "cell_type": "markdown", "id": "a0194e73-dc5e-479e-9940-886f185b9918", "metadata": {}, "source": [ "### Rescaling Temperature Data\n", "\n", "Landsat Collection 2 Level 2 includes measures of [Surface Temperature](https://www.usgs.gov/landsat-missions/landsat-surface-temperature). We'll look at the surface temperature band `TIRS_B10` under the key `lwir11`. The raw values are rescaled, so you should scale and offset the data before interpreting it. Use the metadata in the asset's `raster_bands` to find the scale and offset values:" ] }, { "cell_type": "code", "execution_count": 12, "id": "41612f1d-83b3-443c-b39b-9f91855752e7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'unit': 'kelvin',\n", " 'scale': 0.00341802,\n", " 'nodata': 0,\n", " 'offset': 149.0,\n", " 'data_type': 'uint16',\n", " 'spatial_resolution': 30}" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "band_info = selected_item.assets[\"lwir11\"].extra_fields[\"raster:bands\"][0]\n", "band_info" ] }, { "cell_type": "markdown", "id": "074980a4-b887-43ef-9e7c-ea45585d4aa9", "metadata": {}, "source": [ "To go from the raw values to Kelvin, we apply those values:" ] }, { "cell_type": "code", "execution_count": 13, "id": "b8d2100a-e208-4fff-861b-9da9ea7341af", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'lwir11' (y: 5, x: 5)>\n",
       "array([[297.60525554, 297.56765732, 297.56082128, 297.57449336,\n",
       "        297.58474742],\n",
       "       [297.54031316, 297.46853474, 297.43435454, 297.44802662,\n",
       "        297.4958789 ],\n",
       "       [297.50271494, 297.4275185 , 297.35232206, 297.342068  ,\n",
       "        297.41384642],\n",
       "       [297.49929692, 297.41384642, 297.342068  , 297.32839592,\n",
       "        297.37966622],\n",
       "       [297.52664108, 297.46853474, 297.41726444, 297.41384642,\n",
       "        297.45144464]])\n",
       "Coordinates:\n",
       "  * y            (y) float64 5.289e+06 5.289e+06 5.289e+06 5.288e+06 5.288e+06\n",
       "  * x            (x) float64 5.543e+05 5.544e+05 5.544e+05 5.544e+05 5.544e+05\n",
       "    spatial_ref  int32 32610\n",
       "    time         datetime64[ns] 2021-07-25T18:55:39.475647\n",
       "Attributes:\n",
       "    nodata:   0
" ], "text/plain": [ "\n", "array([[297.60525554, 297.56765732, 297.56082128, 297.57449336,\n", " 297.58474742],\n", " [297.54031316, 297.46853474, 297.43435454, 297.44802662,\n", " 297.4958789 ],\n", " [297.50271494, 297.4275185 , 297.35232206, 297.342068 ,\n", " 297.41384642],\n", " [297.49929692, 297.41384642, 297.342068 , 297.32839592,\n", " 297.37966622],\n", " [297.52664108, 297.46853474, 297.41726444, 297.41384642,\n", " 297.45144464]])\n", "Coordinates:\n", " * y (y) float64 5.289e+06 5.289e+06 5.289e+06 5.288e+06 5.288e+06\n", " * x (x) float64 5.543e+05 5.544e+05 5.544e+05 5.544e+05 5.544e+05\n", " spatial_ref int32 32610\n", " time datetime64[ns] 2021-07-25T18:55:39.475647\n", "Attributes:\n", " nodata: 0" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "temperature = data[\"lwir11\"].astype(float)\n", "temperature *= band_info[\"scale\"]\n", "temperature += band_info[\"offset\"]\n", "temperature[:5, :5]" ] }, { "cell_type": "markdown", "id": "4a48a225-ce61-43e2-b59d-e2f4fea62038", "metadata": {}, "source": [ "To convert from Kelvin to degrees, subtract 273.15." ] }, { "cell_type": "code", "execution_count": 14, "id": "967d43f6-1e76-4822-813e-8bbcfd7639be", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "celsius = temperature - 273.15\n", "celsius.plot(cmap=\"magma\", size=10);" ] }, { "cell_type": "markdown", "id": "1fa0b57e", "metadata": {}, "source": [ "### Landsat Collection 2 Level-1 Data\n", "\n", "Thus far we have worked with Landsat Collection 2 Level-2 data, which is processed to a consistent set of surface reflectance and surface temperature [science products](https://www.usgs.gov/landsat-missions/landsat-science-products). \n", "\n", "The Planetary Computer also hosts Collection 2 Level-1 data (top of atmosphere values) acquired by the [Multispectral Scanner System](https://landsat.gsfc.nasa.gov/multispectral-scanner-system/) (MSS) onboard Landsat 1 through 5. These data do not include a blue band, so a natural color image is not possible. We will plot a color infrared image from the nir08, red, and green bands instead.\n", "\n", "As before, we use pystac-client to search over the `landsat-c2-l1` collection. We'll use the same area of interest and return only those with less than 10% cloud cover." ] }, { "cell_type": "code", "execution_count": 15, "id": "5235cd43", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Returned 272 Items\n" ] } ], "source": [ "search = catalog.search(\n", " collections=[\"landsat-c2-l1\"],\n", " bbox=bbox_of_interest,\n", " query={\"eo:cloud_cover\": {\"lt\": 10}},\n", ")\n", "\n", "items = search.item_collection()\n", "print(f\"Returned {len(items)} Items\")" ] }, { "cell_type": "markdown", "id": "fc09a11a", "metadata": {}, "source": [ "Choose the least cloudy Item." ] }, { "cell_type": "code", "execution_count": 16, "id": "c2995560", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Choosing LM05_L1TP_046027_20121004_02_T2 from 2012-10-04 with 0.0% cloud cover\n" ] } ], "source": [ "selected_item = min(items, key=lambda item: eo.ext(item).cloud_cover)\n", "\n", "print(\n", " f\"Choosing {selected_item.id} from {selected_item.datetime.date()}\"\n", " + f\" with {selected_item.properties['eo:cloud_cover']}% cloud cover\"\n", ")" ] }, { "cell_type": "markdown", "id": "2d9504fd", "metadata": {}, "source": [ "### Available assets\n", "\n", "MSS data has four EO bands assets and a number of metadata files and bands." ] }, { "cell_type": "code", "execution_count": 17, "id": "0df4d6a1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " red: Red Band\n", " green: Green Band\n", " nir08: Near Infrared Band 0.8\n", " nir09: Near Infrared Band 0.9\n", " mtl.txt: Product Metadata File (txt)\n", " mtl.xml: Product Metadata File (xml)\n", " mtl.json: Product Metadata File (json)\n", " qa_pixel: Pixel Quality Assessment Band\n", " qa_radsat: Radiometric Saturation and Dropped Pixel Quality Assessment Band\n", " tilejson: TileJSON with default rendering\n", "rendered_preview: Rendered preview\n" ] } ], "source": [ "max_key_length = len(max(selected_item.assets, key=len))\n", "for key, asset in selected_item.assets.items():\n", " print(f\"{key.rjust(max_key_length)}: {asset.title}\")" ] }, { "cell_type": "markdown", "id": "9a67ee1a", "metadata": {}, "source": [ "We'll use `odc-stac` to load only the EO bands and area we are interested in." ] }, { "cell_type": "code", "execution_count": 18, "id": "1eae9e4d", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "bands_of_interest = [\"nir08\", \"red\", \"green\"]\n", "data = odc.stac.stac_load(\n", " [selected_item], bands=bands_of_interest, bbox=bbox_of_interest\n", ").isel(time=0)\n", "\n", "cir = data.to_array()\n", "\n", "fig, ax = plt.subplots(figsize=(10, 10))\n", "cir.plot.imshow(robust=True, ax=ax)\n", "ax.set_title(\"Color Infrared, Redmond, WA\");" ] } ], "metadata": { "interpreter": { "hash": "ed3fd30edefd633ef3d52b49fba6f05dba1d7f4cc29f008d827ab52268d3b96c" }, "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.10.6" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }