{ "cells": [ { "cell_type": "markdown", "id": "2c0e6328", "metadata": {}, "source": [ "## Accessing ESA WorldCover classification data with the Planetary Computer STAC API\n", "\n", "The European Space Agency (ESA) [WorldCover](https://esa-worldcover.org/en) product provides global land cover maps for the years 2020 and 2021 at 10 meter resolution based on the combination of [Sentinel-1](https://sentinel.esa.int/web/sentinel/missions/sentinel-1) radar data and [Sentinel-2](https://sentinel.esa.int/web/sentinel/missions/sentinel-2) imagery. The discrete classification maps provide 11 classes defined using the Land Cover Classification System (LCCS) developed by the United Nations (UN) Food and Agriculture Organization (FAO). The map images are stored in [cloud-optimized GeoTIFF](https://www.cogeo.org/) format.\n", "\n", "Two versions of the WorldCover product are available:\n", "\n", "- WorldCover 2020 produced using v100 of the algorithm ([User Manual](https://esa-worldcover.s3.eu-central-1.amazonaws.com/v100/2020/docs/WorldCover_PUM_V1.0.pdf))\n", "- WorldCover 2021 produced using v200 of the algorithm ([User Manual](https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/docs/WorldCover_PUM_V2.0.pdf))\n", "\n", "Since the WorldCover maps for 2020 and 2021 were generated with different algorithm versions (v100 and v200, respectively), changes between the maps include both changes in real land cover and changes due to the used algorithms." ] }, { "cell_type": "markdown", "id": "a0fbf056", "metadata": {}, "source": [ "### Data Access\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. If you are using the [Planetary Computer Hub](https://planetarycomputer.microsoft.com/compute) to run this notebook, then your API key is automatically set to the environment variable `PC_SDK_SUBSCRIPTION_KEY` for you when your server is started. Otherwise, you can view your keys by signing in to the [developer portal](https://planetarycomputer.developer.azure-api.net/). The API key may be manually set via the environment variable `PC_SDK_SUBSCRIPTION_KEY` or the following code:\n", "\n", "```python\n", "import planetary_computer\n", "planetary_computer.settings.set_subscription_key()\n", "```" ] }, { "cell_type": "markdown", "id": "36147250", "metadata": {}, "source": [ "The datasets hosted by the Planetary Computer are available in [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": "markdown", "id": "198b54e7", "metadata": {}, "source": [ "### Define the area of interest and search the collection" ] }, { "cell_type": "markdown", "id": "a7c31ef2", "metadata": {}, "source": [ "Let's define a bounding box around Mount Elgon, which sits on the border of Uganda and Kenya." ] }, { "cell_type": "code", "execution_count": 1, "id": "d1a949e6", "metadata": { "tags": [] }, "outputs": [], "source": [ "bbox_of_interest = [33.984, 0.788, 34.902, 1.533]" ] }, { "cell_type": "markdown", "id": "fa58f1c9", "metadata": {}, "source": [ "Use [pystac-client](https://github.com/stac-utils/pystac-client) to search over the ESA WorldCover collection." ] }, { "cell_type": "code", "execution_count": 2, "id": "a56a65b0", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pystac_client\n", "import planetary_computer\n", "\n", "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")\n", "\n", "search = catalog.search(\n", " collections=[\"esa-worldcover\"],\n", " bbox=bbox_of_interest,\n", ")\n", "\n", "items = list(search.get_items())\n", "items" ] }, { "cell_type": "markdown", "id": "b10b12e3", "metadata": {}, "source": [ "### Available Assets & Metadata" ] }, { "cell_type": "markdown", "id": "a5754154", "metadata": {}, "source": [ "Our search returned one item, a 3x3 degree tile of classification data. Let's display the available assets and metadata. " ] }, { "cell_type": "code", "execution_count": 3, "id": "18e5530e", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃ Key              ┃ Value                           ┃\n",
       "┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ map              │ Land Cover Classes              │\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[1mValue \u001b[0m\u001b[1m \u001b[0m┃\n", "┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n", "│ map │ Land Cover Classes │\n", "│ tilejson │ TileJSON with default rendering │\n", "│ rendered_preview │ Rendered preview │\n", "└──────────────────┴─────────────────────────────────┘\n" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import rich.table\n", "\n", "# Assets\n", "t_assets = rich.table.Table(\"Key\", \"Value\")\n", "for key, asset in items[0].assets.items():\n", " t_assets.add_row(key, asset.title)\n", "t_assets" ] }, { "cell_type": "code", "execution_count": 4, "id": "7e6cd3a6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n",
       "┃ Key                         ┃ Value                                                  ┃\n",
       "┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩\n",
       "│ created                     │ 2022-05-16T16:37:31.807816Z                            │\n",
       "│ datetime                    │ None                                                   │\n",
       "│ description                 │ ESA WorldCover product at 10m resolution for year 2020 │\n",
       "│ end_datetime                │ 2020-12-31T23:59:59Z                                   │\n",
       "│ esa_worldcover:product_tile │ N00E033                                                │\n",
       "│ instruments                 │ ['c-sar', 'msi']                                       │\n",
       "│ mission                     │ sentinel-1, sentinel-2                                 │\n",
       "│ platform                    │ sentinel-1a, sentinel-1b, sentinel-2a, sentinel-2b     │\n",
       "│ proj:epsg                   │ 4326                                                   │\n",
       "│ start_datetime              │ 2020-01-01T00:00:00Z                                   │\n",
       "└─────────────────────────────┴────────────────────────────────────────────────────────┘\n",
       "
\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-05-16T16:37:31.807816Z │\n", "│ datetime │ None │\n", "│ description │ ESA WorldCover product at 10m resolution for year 2020 │\n", "│ end_datetime │ 2020-12-31T23:59:59Z │\n", "│ esa_worldcover:product_tile │ N00E033 │\n", "│ instruments │ ['c-sar', 'msi'] │\n", "│ mission │ sentinel-1, sentinel-2 │\n", "│ platform │ sentinel-1a, sentinel-1b, sentinel-2a, sentinel-2b │\n", "│ proj:epsg │ 4326 │\n", "│ start_datetime │ 2020-01-01T00:00:00Z │\n", "└─────────────────────────────┴────────────────────────────────────────────────────────┘\n" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Metadata\n", "t_metadata = rich.table.Table(\"Key\", \"Value\")\n", "for k, v in sorted(items[0].properties.items()):\n", " t_metadata.add_row(k, str(v))\n", "t_metadata" ] }, { "cell_type": "markdown", "id": "99e3ad85", "metadata": {}, "source": [ "The `map` asset contains the full resolution data, but we can use the `rendered_preview` asset to quickly visualize the Item data using the Planetary Computer's data API." ] }, { "cell_type": "code", "execution_count": 5, "id": "65fdf67e", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import Image\n", "\n", "Image(url=items[0].assets[\"rendered_preview\"].href)" ] }, { "cell_type": "markdown", "id": "a3673c18", "metadata": {}, "source": [ "### Render our area of interest\n" ] }, { "cell_type": "markdown", "id": "be28f708", "metadata": {}, "source": [ "Let's create a plot of our area of interest. This dataset includes classification information which maps raster values to class descriptions and preferred colors. We'll want to use this information when plotting the data. We can extract it from the `map` asset." ] }, { "cell_type": "code", "execution_count": 6, "id": "fc56dbc6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┏━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━┓\n",
       "┃ Value ┃ Description              ┃ Hex Color ┃\n",
       "┡━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━┩\n",
       "│ 10    │ Tree cover               │ 006400    │\n",
       "│ 20    │ Shrubland                │ FFBB22    │\n",
       "│ 30    │ Grassland                │ FFFF4C    │\n",
       "│ 40    │ Cropland                 │ F096FF    │\n",
       "│ 50    │ Built-up                 │ FA0000    │\n",
       "│ 60    │ Bare / sparse vegetation │ B4B4B4    │\n",
       "│ 70    │ Snow and ice             │ F0F0F0    │\n",
       "│ 80    │ Permanent water bodies   │ 0064C8    │\n",
       "│ 90    │ Herbaceous wetland       │ 0096A0    │\n",
       "│ 95    │ Mangroves                │ 00CF75    │\n",
       "│ 100   │ Moss and lichen          │ FAE6A0    │\n",
       "└───────┴──────────────────────────┴───────────┘\n",
       "
\n" ], "text/plain": [ "┏━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━┓\n", "┃\u001b[1m \u001b[0m\u001b[1mValue\u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mDescription \u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mHex Color\u001b[0m\u001b[1m \u001b[0m┃\n", "┡━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━┩\n", "│ 10 │ Tree cover │ 006400 │\n", "│ 20 │ Shrubland │ FFBB22 │\n", "│ 30 │ Grassland │ FFFF4C │\n", "│ 40 │ Cropland │ F096FF │\n", "│ 50 │ Built-up │ FA0000 │\n", "│ 60 │ Bare / sparse vegetation │ B4B4B4 │\n", "│ 70 │ Snow and ice │ F0F0F0 │\n", "│ 80 │ Permanent water bodies │ 0064C8 │\n", "│ 90 │ Herbaceous wetland │ 0096A0 │\n", "│ 95 │ Mangroves │ 00CF75 │\n", "│ 100 │ Moss and lichen │ FAE6A0 │\n", "└───────┴──────────────────────────┴───────────┘\n" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "class_list = items[0].assets[\"map\"].extra_fields[\"classification:classes\"]\n", "classmap = {\n", " c[\"value\"]: {\"description\": c[\"description\"], \"hex\": c[\"color-hint\"]}\n", " for c in class_list\n", "}\n", "\n", "t = rich.table.Table(\"Value\", \"Description\", \"Hex Color\")\n", "for k, v in classmap.items():\n", " t.add_row(str(k), v[\"description\"], v[\"hex\"])\n", "t" ] }, { "cell_type": "markdown", "id": "af0af56f", "metadata": {}, "source": [ "We'll use this class information to create a custom colormap for plotting the data and to extract data that we need to add an informative colorbar to the plot." ] }, { "cell_type": "code", "execution_count": 7, "id": "edd12b5b", "metadata": {}, "outputs": [], "source": [ "import matplotlib.colors\n", "\n", "colors = [\"#000000\" for r in range(256)]\n", "for key, value in classmap.items():\n", " colors[int(key)] = f\"#{value['hex']}\"\n", "cmap = matplotlib.colors.ListedColormap(colors)\n", "\n", "# sequences needed for an informative colorbar\n", "values = [key for key in classmap]\n", "boundaries = [(values[i + 1] + values[i]) / 2 for i in range(len(values) - 1)]\n", "boundaries = [0] + boundaries + [255]\n", "ticks = [(boundaries[i + 1] + boundaries[i]) / 2 for i in range(len(boundaries) - 1)]\n", "tick_labels = [value[\"description\"] for value in classmap.values()]" ] }, { "cell_type": "markdown", "id": "50d98429", "metadata": {}, "source": [ "Load the data from our area of interest using `odc-stac`." ] }, { "cell_type": "code", "execution_count": 8, "id": "4462ea97", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'map' (latitude: 7450, longitude: 9180)>\n",
       "array([[30, 30, 30, ..., 30, 30, 30],\n",
       "       [30, 30, 30, ..., 30, 20, 20],\n",
       "       [30, 30, 30, ..., 30, 30, 30],\n",
       "       ...,\n",
       "       [40, 40, 40, ..., 40, 40, 40],\n",
       "       [40, 40, 40, ..., 40, 40, 40],\n",
       "       [40, 40, 40, ..., 40, 40, 40]], dtype=uint8)\n",
       "Coordinates:\n",
       "  * latitude     (latitude) float64 1.533 1.533 1.533 ... 0.7882 0.7881 0.788\n",
       "  * longitude    (longitude) float64 33.98 33.98 33.98 33.98 ... 34.9 34.9 34.9\n",
       "    spatial_ref  int32 4326\n",
       "    time         datetime64[ns] 2020-01-01\n",
       "Attributes:\n",
       "    nodata:   0
" ], "text/plain": [ "\n", "array([[30, 30, 30, ..., 30, 30, 30],\n", " [30, 30, 30, ..., 30, 20, 20],\n", " [30, 30, 30, ..., 30, 30, 30],\n", " ...,\n", " [40, 40, 40, ..., 40, 40, 40],\n", " [40, 40, 40, ..., 40, 40, 40],\n", " [40, 40, 40, ..., 40, 40, 40]], dtype=uint8)\n", "Coordinates:\n", " * latitude (latitude) float64 1.533 1.533 1.533 ... 0.7882 0.7881 0.788\n", " * longitude (longitude) float64 33.98 33.98 33.98 33.98 ... 34.9 34.9 34.9\n", " spatial_ref int32 4326\n", " time datetime64[ns] 2020-01-01\n", "Attributes:\n", " nodata: 0" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import odc.stac\n", "\n", "ds = odc.stac.load(items, crs=\"EPSG:4326\", resolution=0.0001, bbox=bbox_of_interest)\n", "map_data = ds[\"map\"].isel(time=-1).load()\n", "map_data" ] }, { "cell_type": "markdown", "id": "7dbd74df", "metadata": {}, "source": [ "Finally, we can plot our area of interest with a colorbar labeled with the class descriptions." ] }, { "cell_type": "code", "execution_count": 9, "id": "8c201c24-eb03-46ca-8391-b4bd962dbdc0", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from matplotlib import cm\n", "import matplotlib.pyplot as plt\n", "\n", "fig, ax = plt.subplots(figsize=(16, 14))\n", "normalizer = matplotlib.colors.Normalize(vmin=0, vmax=255)\n", "\n", "map_data.isel(latitude=slice(3000, 6000), longitude=slice(4000, 7000)).plot(\n", " ax=ax, cmap=cmap, norm=normalizer\n", ")\n", "\n", "colorbar = fig.colorbar(\n", " cm.ScalarMappable(norm=normalizer, cmap=cmap),\n", " boundaries=boundaries,\n", " values=values,\n", " cax=fig.axes[1].axes,\n", ")\n", "colorbar.set_ticks(ticks, labels=tick_labels)\n", "\n", "ax.set_axis_off()\n", "ax.set_title(\"ESA WorldCover at Mount Elgin\");" ] } ], "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.10.9" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }