{ "cells": [ { "cell_type": "markdown", "id": "4c89ab28-5a6d-427e-8e8e-d9315395e440", "metadata": {}, "source": [ "## Customizable radiometric terrain corrections of Sentinel-1 products\n", "\n", "The Planetary Computer includes both [Sentinel-1 Level-1 Ground Range Detected (GRD)](https://planetarycomputer.microsoft.com/dataset/sentinel-1-grd) and [Sentinel-1 Radiometrically Terrain Corrected (RTC)](https://planetarycomputer.microsoft.com/dataset/sentinel-1-rtc) collections. This tutorial explains the some background on Synthetic Aperture Radar (SAR), geometric and radiometric terrain correction, and introduces [sarsen](https://github.com/bopen/sarsen), an open-source library for working with SAR data. `sarsen` enables geometric and radiometric corrections using any digital elevation model (DEM) supported by GDAL / Proj.\n", "\n", "### Background\n", "\n", "The typical side-looking [Synthetic Aperture Radar (SAR)](https://earthdata.nasa.gov/learn/backgrounders/what-is-sar) system acquires data with uniform sampling in azimuth and slant range, where the azimuth and range represents the time when a given target is acquired and the absolute sensor-to-target distance, respectively.\n", "Because of this, the near range appears compressed with respect to the far range. Furthermore, any deviation of the target elevation from a smooth geoid results in additional local geometric and radiometric distortions known as foreshortening, layover and shadow.\n", "\n", "- Radar foreshortening: Terrain surfaces sloping towards the radar appear shortened relative to those sloping away from the radar.\n", "- Radar layover: It's an extreme case of foreshortening occurring when the terrain slope is greater than the angle of the incident signal. \n", "- Radar shadows: They occur when ground points at the same azimuth but different slant ranges are aligned in the direction of the line-of-sight. This is usually due to a back slope with an angle steeper than the viewing angle. When this happens, the radar signal never reaches the farthest points, and thus there is no measurement, meaning that this lack of information is unrecoverable.\n", "\n", "### Geometric Terrain Correction\n", "\n", "The [Sentinel-1 GRD](https://planetarycomputer.microsoft.com/dataset/sentinel-1-grd) product already provides a geometric correction that removes the compression effect on the near-range. Geometric Terrain Correction *also* corrects for displacements due to target elevation. When applying geometric terrain correction, the resolution and accuracy of the input Digital Elevation Model (DEM). The Planetary Computer has a number of [DEMs](https://planetarycomputer.microsoft.com/catalog?tags=DEM) to choose from.\n", "\n", "### Radiometric Terrain Correction\n", "\n", "Terrain variations affect both the position of a given point on the Earth's surface *and* the brightness of the radar return. The `sarsen` radiometric terrain correction compensates for the backscatter modulation generated by the topography of the scene. This produces a more uniform backscatter image, emphasizing the radiometric differences of the terrain.\n", "\n", "The accuracy of Radiometric Terrain Corrected (RTC) products is also strongly affected by the resolution of the input DEM.\n", "\n", "This tutorial shows how to perform (i) geometric and (ii) radiometric terrain corrections on the Sentinel-1 GRD product using `sarsen`. We use a 10-meter resolution DEM, the same resolution of the DEM used to generate the [Sentinel-1 RTC product](https://planetarycomputer.microsoft.com/dataset/sentinel-1-rtc) available on the Planetary Computer. The comparison at the end of this notebook demonstrates that the RTC computed by `sarsen` is consistent with the RTC from the Planetary Computer.\n", "\n", "As an example, we use data covering the South-of-Redmond region (Seattle, US).\n", "\n", "Steps:\n", "\n", "- Download the Sentinel-1 GRD\n", "- Download the 10-meter DEM\n", "- Compute the GTC using `sarsen`\n", "- Compute the RTC using `sarsen`\n", "- Compare the GTC to the RTC\n", "- Compare the RTC computed using `sarsen` to the RTC already available on the Planetery Computer \n", "\n", "**Note**: Download/retrieval steps are slower on local machines compared to the Planetary Computer. In future versions, it will be possible to access data via [fsspec](https://filesystem-spec.readthedocs.io/en/latest/) without having to download data locally." ] }, { "cell_type": "code", "execution_count": 1, "id": "f95d095c-a5b6-40d2-a47c-1adc4e60d31b", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "plt.rcParams[\"figure.figsize\"] = (10, 7)" ] }, { "cell_type": "code", "execution_count": 2, "id": "c949e993-be36-42c6-87f0-773093a464a0", "metadata": {}, "outputs": [], "source": [ "import os\n", "import tempfile\n", "\n", "import rioxarray\n", "\n", "from sarsen import apps\n", "\n", "import adlfs\n", "import planetary_computer\n", "import pystac_client\n", "import stackstac" ] }, { "cell_type": "markdown", "id": "23a700d1-2b4a-45ba-8a64-bd154f775b76", "metadata": {}, "source": [ "### Processing definitions" ] }, { "cell_type": "code", "execution_count": 3, "id": "2cfc92e0-91ee-47f6-8329-f002eb6409b1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'/tmp'" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# create a temporary directory where to store downloaded data\n", "tmp_dir = tempfile.gettempdir()\n", "\n", "# DEM path\n", "dem_path = os.path.join(tmp_dir, \"South-of-Redmond-10m.tif\")\n", "\n", "# path to Sentinel-1 input product in the Planetary Computer\n", "product_folder = \"GRD/2021/12/17/IW/DV/S1B_IW_GRDH_1SDV_20211217T141304_20211217T141329_030066_039705_9048\" # noqa: E501\n", "\n", "# band to be processed\n", "measurement_group = \"IW/VV\"\n", "\n", "tmp_dir" ] }, { "cell_type": "markdown", "id": "868bd066-4122-4f5c-932e-8b607388ae1f", "metadata": {}, "source": [ "#### Area of interest definition: South-of-Redmond (Seattle, US)" ] }, { "cell_type": "code", "execution_count": 4, "id": "51cfd0f3-1be4-493a-a419-b56425c139eb", "metadata": {}, "outputs": [], "source": [ "lon, lat = [-121.95, 47.04]\n", "buffer = 0.2\n", "bbox = [lon - buffer, lat - buffer, lon + buffer, lat + buffer]" ] }, { "cell_type": "markdown", "id": "211d1633-cd35-4e1c-8b71-bc6432e7495c", "metadata": {}, "source": [ "#### DEMs discovery\n", "\n", "Here we use the [USGS 3dep-seamles DEM](https://planetarycomputer.microsoft.com/dataset/3dep-seamless) with a 10-meter ground sample distance (GDS). Note that **any DEM supported by GDAL/Proj can be used**.\n", "\n", "Using `pystac_client` we can search the Planetary Computer's STAC endpoint for items matching our query parameters. \n", "As multiple DEMs acquired at different times are available in this area, we select the DEMs with 10-meter GDS and perform the average of the remaining DEMs along the time dimension." ] }, { "cell_type": "code", "execution_count": 5, "id": "8a8883ba-e53f-41bf-95a3-25219c1799e6", "metadata": {}, "outputs": [], "source": [ "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")\n", "search = catalog.search(\n", " collections=\"3dep-seamless\", bbox=bbox, query={\"gsd\": {\"eq\": 10}}\n", ")\n", "items = search.item_collection()" ] }, { "cell_type": "markdown", "id": "0d016965", "metadata": {}, "source": [ "Here we load the data into an xarray `DataArray` using stackstac." ] }, { "cell_type": "code", "execution_count": 6, "id": "e7021623-ba3f-407e-8c96-f09824b2adbe", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray 'stackstac-91c9486d3e17af5ae0ea1ccefd39911d' (time: 4,\n", " y: 4321, x: 4321)>\n", "dask.array<getitem, shape=(4, 4321, 4321), dtype=float64, chunksize=(1, 1024, 1024), chunktype=numpy.ndarray>\n", "Coordinates: (12/13)\n", " * time (time) datetime64[ns] 2018-02-02 2018-02-08 ... 2020-01-07\n", " id (time) <U10 'n48w122-13' 'n47w122-13' ... 'n48w123-13'\n", " band <U4 'data'\n", " * x (x) float64 -122.2 -122.1 -122.1 ... -121.8 -121.8 -121.8\n", " * y (y) float64 47.24 47.24 47.24 47.24 ... 46.84 46.84 46.84\n", " threedep:region <U7 'n40w130'\n", " ... ...\n", " proj:epsg int64 5498\n", " proj:shape object {10812}\n", " gsd int64 10\n", " end_datetime (time) <U20 '2016-12-31T00:00:00Z' ... '2019-04-25T00:00...\n", " description <U1849 'This tile of the 3D Elevation Program (3DEP) sea...\n", " epsg int64 5498\n", "Attributes:\n", " spec: RasterSpec(epsg=5498, bounds=(-122.15000053746, 46.839907613...\n", " crs: epsg:5498\n", " transform: | 0.00, 0.00,-122.15|\\n| 0.00,-0.00, 47.24|\\n| 0.00, 0.00, 1...\n", " resolution: 9.2592593e-05
<xarray.DataArray 'stackstac-d90c33d750479ce84c6f5eda1b1f3eec' (y: 3000, x: 2900)>\n", "array([[ 196.68109131, 196.14202881, 194.11082458, ..., 820.71972656,\n", " 820.05755615, 819.68597412],\n", " [ 196.68109131, 196.14202881, 195.13980103, ..., 817.84838867,\n", " 814.0760498 , 814.61773682],\n", " [ 197.03511047, 195.65205383, 195.13980103, ..., 812.12884521,\n", " 810.13116455, 809.48669434],\n", " ...,\n", " [ 602.28057861, 607.78918457, 609.63909912, ..., 4130.93261719,\n", " 4135.22167969, 4137.34033203],\n", " [ 597.96710205, 602.36712646, 603.93676758, ..., 4133.00195312,\n", " 4137.37695312, 4139.43896484],\n", " [ 594.58496094, 599.09234619, 600.71466064, ..., 4135.22558594,\n", " 4139.52978516, 4141.82373047]])\n", "Coordinates:\n", " * x (x) float64 5.65e+05 5.65e+05 ... 5.94e+05 5.94e+05\n", " * y (y) float64 5.22e+06 5.22e+06 ... 5.19e+06 5.19e+06\n", " proj:shape object {10812}\n", " band <U4 'data'\n", " threedep:region <U7 'n40w130'\n", " epsg int64 5498\n", " proj:epsg int64 5498\n", " description <U1849 'This tile of the 3D Elevation Program (3DEP) sea...\n", " gsd int64 10\n", " spatial_ref int64 0\n", "Attributes:\n", " _FillValue: 1.7976931348623157e+308" ], "text/plain": [ "