{ "cells": [ { "cell_type": "markdown", "id": "4ee7f9a2", "metadata": {}, "source": [ "## Coregistration of Satellite Image Data Sources\n", "\n", "In this tutorial, you'll learn how to use [stackstac](https://stackstac.readthedocs.io/en/latest/) on Planetary Computer to coregister (align) two spatially overlapping raster images from Sentinel-2 Level 2A and Landsat Collection 2 Level-2 into a single dataset. You will then calculate the [Normalized difference vegetation index (NDVI)](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index) of the resulting dataset and compare it to the NDVI based on the original data.\n", "\n", "This tutorial covers the following steps:\n", "\n", "- Searching for and loading data from the Planetary Computer's [STAC endpoint](../quickstarts/reading-stac.ipynb)\n", "- Reprojecting and resampling [Sentinel-2 Level-2A](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a) and [Landsat Collection 2 Level-2\n", "](https://planetarycomputer.microsoft.com/dataset/landsat-c2-l2) data to use a common resolution and coordinate system\n", "- Croping to the region of interest and aligning Sentinel and Landsat data into a single dataset\n", "- Calculating and comparing NDVIs with [xarray-spatial](https://xarray-spatial.org/index.html)" ] }, { "cell_type": "code", "execution_count": 1, "id": "80ea0978", "metadata": {}, "outputs": [], "source": [ "import planetary_computer\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "import stackstac\n", "import pystac_client\n", "import rasterio\n", "import xrspatial.multispectral as ms" ] }, { "cell_type": "markdown", "id": "b4ace20e", "metadata": {}, "source": [ "### Preparation: create a local Dask cluster\n", "\n", "In this tutorial, you'll be using a small dataset. Create a local Dask cluster to hold your data." ] }, { "cell_type": "code", "execution_count": 2, "id": "d5a3a208-184a-464f-ae28-533ef03e460d", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-08-16 16:22:54,567 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-fh4atg6f', purging\n", "2022-08-16 16:22:54,567 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-q6v1o3c3', purging\n", "2022-08-16 16:22:54,567 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-ces756c5', purging\n", "2022-08-16 16:22:54,567 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-beygsd3u', purging\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "/proxy/8787/status\n" ] } ], "source": [ "from dask.distributed import Client\n", "\n", "client = Client()\n", "print(f\"/proxy/{client.scheduler_info()['services']['dashboard']}/status\")" ] }, { "cell_type": "markdown", "id": "2d47e6bc", "metadata": {}, "source": [ "To follow the progress of your computations, you can [access the Dask Dashboard](https://planetarycomputer.microsoft.com/docs/quickstarts/scale-with-dask/#Open-the-dashboard) at the URL from the previous cell's output.\n", "\n", "### Load matching Sentinel and Landsat datasets\n", "\n", "The area of interest covers a small section of the Mississippi River, USA, near the Yazoo National Wildlife Refuge. The two datasets have different tiling boundaries, so scenes from one won't line up perfectly with scenes from the other.\n", "\n", "We'll select a single Sentinel scene and then make a median mosaic of Landsat scenes overlapping with that Sentinel scene." ] }, { "cell_type": "code", "execution_count": 3, "id": "397ee012", "metadata": {}, "outputs": [], "source": [ "area_of_interest = {\n", " \"type\": \"Polygon\",\n", " \"coordinates\": [\n", " [\n", " [-91.92, 33.43],\n", " [-90.74, 33.41],\n", " [-90.76, 32.42],\n", " [-91.93, 32.44],\n", " [-91.92, 33.43],\n", " ]\n", " ],\n", "}\n", "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")\n", "\n", "sentinel_search = catalog.search(\n", " intersects=area_of_interest,\n", " datetime=\"2020-09-07/2020-09-08\",\n", " collections=[\"sentinel-2-l2a\"],\n", ")\n", "\n", "sentinel_item = next(sentinel_search.items()) # select the first item" ] }, { "cell_type": "markdown", "id": "ab0a46ed", "metadata": {}, "source": [ "Now let's, use [pystac_client](https://pystac-client.readthedocs.io/en/latest/) to identify all matching Landsat images with a cloud coverage of less than 10 percent in September 2020 that overlap with that Sentinel scene's bounding box." ] }, { "cell_type": "code", "execution_count": 4, "id": "d03bd2d8", "metadata": {}, "outputs": [], "source": [ "time_of_interest = \"2020-09-01/2020-09-30\"\n", "\n", "search = catalog.search(\n", " collections=[\"landsat-c2-l2\"],\n", " intersects=sentinel_item.geometry,\n", " datetime=time_of_interest,\n", " query={\n", " \"eo:cloud_cover\": {\n", " \"lt\": 10,\n", " },\n", " \"platform\": {\"in\": [\"landsat-8\", \"landsat-9\"]},\n", " },\n", ")\n", "\n", "landsat_items = search.item_collection()" ] }, { "cell_type": "markdown", "id": "07d1f6dd", "metadata": {}, "source": [ "### Match Sentinel and Landsat data properties\n", "\n", "Let's start by loading the Sentinel data into an `xarray.DataArray` using its native [Coordinate Reference System (CRS)](https://gdal.org/tutorials/osr_api_tut.html) `epsg=32615` and a resolution of 100m." ] }, { "cell_type": "code", "execution_count": 5, "id": "d1c7ec4a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(RasterSpec(epsg=32615, bounds=(699900, 3690200, 809800, 3800100), resolutions_xy=(100, 100)),\n", " 100,\n", " (4, 1099, 1099))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentinel_data = (\n", " (\n", " stackstac.stack(\n", " sentinel_item,\n", " resolution=100,\n", " assets=[\"B02\", \"B03\", \"B04\", \"B08\"], # blue, green, red, nir\n", " )\n", " .where(lambda x: x > 0, other=np.nan) # sentinel-2 uses 0 as nodata\n", " .assign_coords(band=lambda x: x.common_name.rename(\"band\")) # use common names\n", " )\n", " .isel(time=0)\n", " .persist()\n", ")\n", "\n", "sentinel_data.spec, sentinel_data.resolution, sentinel_data.shape" ] }, { "cell_type": "markdown", "id": "4be7bdd1", "metadata": {}, "source": [ "Before coregistering data, you need to make sure that both data sets use the same resolution and the same CRS. Use `stackstac.stack` to set the `epsg` and `resolution` properties of your Landast data to match the CRS and resolution of your Sentinel data as you read it into xarray.\n", "\n", "Translating data from one CRS to another is called reprojection. See [Reprojecting](reprojection.ipynb) to learn more." ] }, { "cell_type": "code", "execution_count": 6, "id": "ab8f06e4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(RasterSpec(epsg=32615, bounds=(699900, 3690200, 809800, 3800100), resolutions_xy=(100, 100)),\n", " (4, 1099, 1099))" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "landsat_data = (\n", " (\n", " stackstac.stack(\n", " landsat_items,\n", " resolution=sentinel_data.resolution, # resample to Sentinel data resolution\n", " epsg=sentinel_data.spec.epsg, # reporoject to CRS of Sentinel data\n", " bounds=sentinel_data.spec.bounds, # set bounds to match Sentinel data\n", " assets=[\"blue\", \"green\", \"red\", \"nir08\"],\n", " resampling=rasterio.enums.Resampling.bilinear,\n", " ).where(\n", " lambda x: x > 0, other=np.nan\n", " ) # landsat-c2-l2 uses 0 as nodata\n", " )\n", " .median(dim=\"time\", keep_attrs=True)\n", " .persist()\n", ")\n", "\n", "landsat_data.spec, landsat_data.shape" ] }, { "cell_type": "markdown", "id": "8c8cc752", "metadata": {}, "source": [ "Both the Sentinel and the Landsat data should now have the same raster spec information. To make sure, check whether their x and y coordinates match." ] }, { "cell_type": "code", "execution_count": 7, "id": "f8f3395b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(landsat_data.x.data == sentinel_data.x.data).all()" ] }, { "cell_type": "code", "execution_count": 8, "id": "bff09a1d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(landsat_data.y.data == sentinel_data.y.data).all()" ] }, { "cell_type": "markdown", "id": "226cfa5d", "metadata": {}, "source": [ "To display a preview of the two data sets, use the [xrspatial.multispectral.true_color](https://xarray-spatial.org/reference/_autosummary/xrspatial.multispectral.true_color.html) function with `sentinel_data` and `landsat_data`." ] }, { "cell_type": "code", "execution_count": 9, "id": "dedf9ab3", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sentinel_image = ms.true_color(\n", " sentinel_data.sel(band=\"red\"),\n", " sentinel_data.sel(band=\"green\"),\n", " sentinel_data.sel(band=\"blue\"),\n", " c=25,\n", " th=0.1,\n", " name=\"Sentinel-2\",\n", ")\n", "\n", "landsat_image = ms.true_color(\n", " landsat_data.sel(band=\"red\"),\n", " landsat_data.sel(band=\"green\"),\n", " landsat_data.sel(band=\"blue\"),\n", " c=15,\n", " th=0.1,\n", " name=\"Landsat C2-L2\",\n", ")\n", "\n", "org_imgs = xr.concat(\n", " [sentinel_image, landsat_image],\n", " pd.Index([sentinel_image.name, landsat_image.name], name=\"name\"),\n", ")\n", "\n", "org_imgs.plot.imshow(x=\"x\", y=\"y\", col=\"name\", figsize=(15, 8));" ] }, { "cell_type": "markdown", "id": "3a9ea845", "metadata": {}, "source": [ "### Coregister Sentinel and Landsat data into one dataset\n", "\n", "Combine Sentinel and Landsat data into a single [xarray Dataset](http://xarray.pydata.org/en/stable/generated/xarray.Dataset.html)." ] }, { "cell_type": "code", "execution_count": 10, "id": "6b47b611", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:   (band: 4, y: 1099, x: 1099)\n",
       "Coordinates:\n",
       "  * x         (x) float64 6.999e+05 7e+05 7.001e+05 ... 8.096e+05 8.097e+05\n",
       "  * y         (y) float64 3.8e+06 3.8e+06 3.8e+06 ... 3.69e+06 3.69e+06\n",
       "  * band      (band) <U5 'blue' 'green' 'red' 'nir'\n",
       "Data variables:\n",
       "    sentinel  (band, y, x) float64 dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>\n",
       "    landsat   (band, y, x) float64 dask.array<chunksize=(4, 1099, 1099), meta=np.ndarray>\n",
       "Attributes:\n",
       "    spec:        RasterSpec(epsg=32615, bounds=(699900, 3690200, 809800, 3800...\n",
       "    crs:         epsg:32615\n",
       "    transform:   | 100.00, 0.00, 699900.00|\\n| 0.00,-100.00, 3800100.00|\\n| 0...\n",
       "    resolution:  100
" ], "text/plain": [ "\n", "Dimensions: (band: 4, y: 1099, x: 1099)\n", "Coordinates:\n", " * x (x) float64 6.999e+05 7e+05 7.001e+05 ... 8.096e+05 8.097e+05\n", " * y (y) float64 3.8e+06 3.8e+06 3.8e+06 ... 3.69e+06 3.69e+06\n", " * band (band) \n", " landsat (band, y, x) float64 dask.array\n", "Attributes:\n", " spec: RasterSpec(epsg=32615, bounds=(699900, 3690200, 809800, 3800...\n", " crs: epsg:32615\n", " transform: | 100.00, 0.00, 699900.00|\\n| 0.00,-100.00, 3800100.00|\\n| 0...\n", " resolution: 100" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coregistration_data = xr.Dataset(\n", " data_vars=dict(\n", " sentinel=([\"band\", \"y\", \"x\"], sentinel_data.data),\n", " landsat=([\"band\", \"y\", \"x\"], landsat_data.data),\n", " ),\n", " coords=dict(\n", " x=sentinel_data.x.data,\n", " y=sentinel_data.y.data,\n", " band=sentinel_data.band.data,\n", " ),\n", " attrs=sentinel_data.attrs, # both dataset have the same attributes\n", ")\n", "\n", "coregistration_data" ] }, { "cell_type": "markdown", "id": "03a3b414", "metadata": {}, "source": [ "### Wrap-up: compare NDVI\n", "\n", "Now let's compute NDVI on the coregistered dataset using [xrspatial.multispectral.ndvi](https://xarray-spatial.org/reference/_autosummary/xrspatial.multispectral.ndvi.html). We'll use the near-infrared band from Sentinel and the red band from Landsat. Because the original datasets both contain Red and near-infrared bands, we can compare the NDVI computed on the co-registered datasets to the NDVIs computed on the original datasets.\n", "\n", "First, compute the NVDIs for `sentinel_ndvi` and `landsat_ndvi`." ] }, { "cell_type": "code", "execution_count": 11, "id": "257de64c", "metadata": {}, "outputs": [], "source": [ "sentinel_ndvi = ms.ndvi(\n", " nir_agg=coregistration_data.sentinel.sel(band=\"nir\"),\n", " red_agg=coregistration_data.sentinel.sel(band=\"red\"),\n", " name=\"sentinel_ndvi\",\n", ")\n", "\n", "landsat_ndvi = ms.ndvi(\n", " nir_agg=coregistration_data.landsat.sel(band=\"nir\"),\n", " red_agg=coregistration_data.landsat.sel(band=\"red\"),\n", " name=\"landsat_ndvi\",\n", ")" ] }, { "cell_type": "markdown", "id": "f07798c6", "metadata": {}, "source": [ "And now the NDVI for the combined dataset, using the near-infrared band from Sentinel-2 and the red band from Landsat C2-L2." ] }, { "cell_type": "code", "execution_count": 12, "id": "dc0c0a09", "metadata": {}, "outputs": [], "source": [ "coregistration_ndvi = ms.ndvi(\n", " nir_agg=coregistration_data.sentinel.sel(band=\"nir\"),\n", " red_agg=coregistration_data.landsat.sel(band=\"red\"),\n", " name=\"coregistration_ndvi\",\n", ")" ] }, { "cell_type": "markdown", "id": "7d780cbf", "metadata": {}, "source": [ "Finally, use xarray's [concat](http://xarray.pydata.org/en/stable/generated/xarray.concat.html) function to combine the three datasets into one array and display the results." ] }, { "cell_type": "code", "execution_count": 13, "id": "ed330f22", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ndvi_imgs = xr.concat(\n", " [sentinel_ndvi, landsat_ndvi, coregistration_ndvi],\n", " pd.Index(\n", " [sentinel_ndvi.name, landsat_ndvi.name, coregistration_ndvi.name], name=\"name\"\n", " ),\n", ")\n", "\n", "ndvi_imgs.name = \"ndvi_images\"\n", "\n", "ndvi_imgs.plot.imshow(x=\"x\", y=\"y\", col=\"name\", figsize=(18, 5), cmap=\"viridis\");" ] } ], "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": { "0a7033878b174ca29ff3a8676ff877ff": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "VBoxModel", "state": { "layout": "IPY_MODEL_6a3963dcebfb498c8e8e624b64486fda" } }, "6a3963dcebfb498c8e8e624b64486fda": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": {} } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }