{ "cells": [ { "cell_type": "markdown", "id": "db8d71aa", "metadata": {}, "source": [ "## Radiometric terrain correction qualitative assessment \n", "\n", "[sarsen](https://github.com/bopen/sarsen) radiometric correction implements the gamma flatting algorithm proposed by Small, David (2011). [Flattening Gamma: Radiometric Terrain Correction for SAR Imagery. Geoscience and Remote Sensing, IEEE Transactions on. 49. 3081 - 3093. 10.1109/TGRS.2011.2120616](https://www.researchgate.net/publication/224230688_Flattening_Gamma_Radiometric_Terrain_Correction_for_SAR_Imagery)\n", "\n", "In this notebook, we perform the same qualitative test done by the authors to evaluate the goodness of the gamma flattening correction. For an introduction to terrain correction, see [Customizable radiometric terrain corrections of Sentinel-1 products\n", "](https://planetarycomputer.microsoft.com/docs/tutorials/customizable-rtc-sentinel1/)." ] }, { "cell_type": "markdown", "id": "6d604a80-9222-46f2-97cb-afe5bfb9bdbd", "metadata": {}, "source": [ "### Install Dependencies and Imports\n", "Additional dependencies: `sarsen==0.9.2`" ] }, { "cell_type": "code", "execution_count": 1, "id": "ce375bf1-9ec8-4aa5-b3ce-cbfe5995c31e", "metadata": {}, "outputs": [], "source": [ "# !pip install sarsen==0.9.2" ] }, { "cell_type": "code", "execution_count": 2, "id": "9f82185a-60a4-4651-8bdf-6fb1ee82c884", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'retina'\n", "\n", "\n", "import matplotlib.pyplot as plt\n", "from matplotlib.colors import Normalize\n", "from matplotlib import cm\n", "\n", "plt.rcParams[\"figure.figsize\"] = (10, 7)\n", "plt.rcParams[\"font.size\"] = 14" ] }, { "cell_type": "code", "execution_count": 3, "id": "1391c47f-d76e-4840-8eed-d547b4784331", "metadata": {}, "outputs": [], "source": [ "import os\n", "import tempfile\n", "\n", "import numpy as np\n", "import xarray as xr\n", "from scipy.stats import gaussian_kde as kde\n", "from sarsen import apps, geocoding, orbit, scene\n", "\n", "import adlfs\n", "import planetary_computer\n", "import pystac_client\n", "import stackstac\n", "import rioxarray # noqa: F401\n", "\n", "plt.rcParams[\"figure.figsize\"] = (10, 7)" ] }, { "cell_type": "markdown", "id": "86970ee8", "metadata": {}, "source": [ "### Processing definitions" ] }, { "cell_type": "code", "execution_count": 4, "id": "107608ba-51f5-4a3d-ab70-bc9fb34eafed", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'/tmp'" ] }, "execution_count": 4, "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": "4dd0a0f2-4476-48fa-bceb-6f5a872c59de", "metadata": {}, "source": [ "#### Area of interest definition: South-of-Redmond (Seattle, US)" ] }, { "cell_type": "code", "execution_count": 5, "id": "8ec8652a-7736-40cc-b8d9-1134493e5888", "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": "b144d7a4-839b-4e41-95ca-6d20418432d9", "metadata": {}, "source": [ "#### DEMs discovery\n", "\n", "Here we use the DEM with a 10-meter ground sample distance (GDS) available on the Planetary Computer. 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": 6, "id": "e08b28b1-a706-4005-b1e8-633a748e7f69", "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": "7503fa77-ba5a-492c-9662-ae2f9b2ba5f4", "metadata": {}, "source": [ "Here we load the data into an xarray `DataArray` using stackstac." ] }, { "cell_type": "code", "execution_count": 7, "id": "625eb3d3-90b0-429b-addc-e0ac59e60942", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'stackstac-6fbfb106cd445ae5ea0f677ec80c3c41' (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",
       "    start_datetime   (time) <U20 '1952-01-01T00:00:00Z' ... '2017-09-14T00:00...\n",
       "    ...               ...\n",
       "    proj:shape       object {10812}\n",
       "    proj:epsg        int64 5498\n",
       "    threedep:region  <U7 'n40w130'\n",
       "    gsd              int64 10\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
" ], "text/plain": [ "\n", "dask.array\n", "Coordinates: (12/13)\n", " * time (time) datetime64[ns] 2018-02-02 2018-02-08 ... 2020-01-07\n", " id (time) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'stackstac-6fbfb106cd445ae5ea0f677ec80c3c41' (y: 4487, x: 3100)>\n",
       "array([[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n",
       "        1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n",
       "       [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n",
       "        1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n",
       "       [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n",
       "        1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n",
       "       ...,\n",
       "       [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n",
       "        1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n",
       "       [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n",
       "        1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n",
       "       [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n",
       "        1.79769313e+308, 1.79769313e+308, 1.79769313e+308]])\n",
       "Coordinates:\n",
       "  * x                (x) float64 5.643e+05 5.643e+05 ... 5.953e+05 5.953e+05\n",
       "  * y                (y) float64 5.233e+06 5.233e+06 ... 5.188e+06 5.188e+06\n",
       "    proj:shape       object {10812}\n",
       "    epsg             int64 5498\n",
       "    threedep:region  <U7 'n40w130'\n",
       "    proj:epsg        int64 5498\n",
       "    gsd              int64 10\n",
       "    band             <U4 'data'\n",
       "    description      <U1849 'This tile of the 3D Elevation Program (3DEP) sea...\n",
       "    spatial_ref      int64 0\n",
       "Attributes:\n",
       "    _FillValue:  1.7976931348623157e+308
" ], "text/plain": [ "\n", "array([[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n", " 1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n", " [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n", " 1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n", " [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n", " 1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n", " ...,\n", " [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n", " 1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n", " [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n", " 1.79769313e+308, 1.79769313e+308, 1.79769313e+308],\n", " [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ...,\n", " 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]])\n", "Coordinates:\n", " * x (x) float64 5.643e+05 5.643e+05 ... 5.953e+05 5.953e+05\n", " * y (y) float64 5.233e+06 5.233e+06 ... 5.188e+06 5.188e+06\n", " proj:shape object {10812}\n", " epsg int64 5498\n", " threedep:region \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'stackstac-6fbfb106cd445ae5ea0f677ec80c3c41' (y: 1500, x: 1400)>\n",
       "array([[ 840.86486816,  837.69940186,  836.6675415 , ..., 1222.60290527,\n",
       "        1224.68151855, 1229.60192871],\n",
       "       [ 847.59912109,  843.43511963,  841.19110107, ..., 1226.77722168,\n",
       "        1229.88330078, 1235.05480957],\n",
       "       [ 853.09313965,  849.24804688,  847.56011963, ..., 1231.57214355,\n",
       "        1235.27197266, 1240.26989746],\n",
       "       ...,\n",
       "       [1057.71154785, 1056.44726562, 1053.34716797, ..., 1967.21203613,\n",
       "        1966.72473145, 1967.17895508],\n",
       "       [1053.78833008, 1052.40124512, 1049.59545898, ..., 1970.66943359,\n",
       "        1969.1439209 , 1969.60314941],\n",
       "       [1049.06298828, 1048.03967285, 1045.55126953, ..., 1973.11730957,\n",
       "        1971.60437012, 1971.93334961]])\n",
       "Coordinates:\n",
       "  * x                (x) float64 5.79e+05 5.79e+05 ... 5.93e+05 5.93e+05\n",
       "  * y                (y) float64 5.21e+06 5.21e+06 ... 5.195e+06 5.195e+06\n",
       "    proj:shape       object {10812}\n",
       "    epsg             int64 5498\n",
       "    threedep:region  <U7 'n40w130'\n",
       "    proj:epsg        int64 5498\n",
       "    gsd              int64 10\n",
       "    band             <U4 'data'\n",
       "    description      <U1849 'This tile of the 3D Elevation Program (3DEP) sea...\n",
       "    spatial_ref      int64 0\n",
       "Attributes:\n",
       "    _FillValue:  1.7976931348623157e+308
" ], "text/plain": [ "\n", "array([[ 840.86486816, 837.69940186, 836.6675415 , ..., 1222.60290527,\n", " 1224.68151855, 1229.60192871],\n", " [ 847.59912109, 843.43511963, 841.19110107, ..., 1226.77722168,\n", " 1229.88330078, 1235.05480957],\n", " [ 853.09313965, 849.24804688, 847.56011963, ..., 1231.57214355,\n", " 1235.27197266, 1240.26989746],\n", " ...,\n", " [1057.71154785, 1056.44726562, 1053.34716797, ..., 1967.21203613,\n", " 1966.72473145, 1967.17895508],\n", " [1053.78833008, 1052.40124512, 1049.59545898, ..., 1970.66943359,\n", " 1969.1439209 , 1969.60314941],\n", " [1049.06298828, 1048.03967285, 1045.55126953, ..., 1973.11730957,\n", " 1971.60437012, 1971.93334961]])\n", "Coordinates:\n", " * x (x) float64 5.79e+05 5.79e+05 ... 5.93e+05 5.93e+05\n", " * y (y) float64 5.21e+06 5.21e+06 ... 5.195e+06 5.195e+06\n", " proj:shape object {10812}\n", " epsg int64 5498\n", " threedep:region " ], "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 469, "width": 666 } }, "output_type": "display_data" } ], "source": [ "dem_raster.plot()\n", "_ = plt.title(\"DEM in UTM coordinates\")" ] }, { "cell_type": "markdown", "id": "4445aa2a-abf4-4162-ab0b-b22ae069e4e4", "metadata": {}, "source": [ "### Define GRD parameters" ] }, { "cell_type": "code", "execution_count": 12, "id": "0c35e4ca", "metadata": {}, "outputs": [], "source": [ "grd_account_name = \"sentinel1euwest\"\n", "grd_storage_container = \"s1-grd\"\n", "grd_product_folder = f\"{grd_storage_container}/{product_folder}\"\n", "grd_local_path = os.path.join(tmp_dir, product_folder)" ] }, { "cell_type": "markdown", "id": "ac695ed8-10bb-4efe-a43e-372c289605e3", "metadata": {}, "source": [ "### Retrieve Sentinel-1 GRD" ] }, { "cell_type": "code", "execution_count": 13, "id": "81012b36", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['s1-grd/GRD/2021/12/17/IW/DV/S1B_IW_GRDH_1SDV_20211217T141304_20211217T141329_030066_039705_9048/manifest.safe']" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "grd_token = planetary_computer.sas.get_token(\n", " grd_account_name, grd_storage_container\n", ").token\n", "grd_fs = adlfs.AzureBlobFileSystem(grd_account_name, credential=grd_token)\n", "grd_fs.ls(f\"{grd_product_folder}/manifest.safe\")" ] }, { "cell_type": "code", "execution_count": 14, "id": "0b1daaa9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/tmp/GRD/2021/12/17/IW/DV/S1B_IW_GRDH_1SDV_20211217T141304_20211217T141329_030066_039705_9048\n" ] } ], "source": [ "grd_fs.get(grd_product_folder, grd_local_path, recursive=True)\n", "!ls -d {grd_local_path}" ] }, { "cell_type": "markdown", "id": "5dbe585a-5467-445e-a9cb-cd8c5b66fe8e", "metadata": {}, "source": [ "### Process" ] }, { "cell_type": "markdown", "id": "54a3d505-a01c-4a88-a939-0169264609d0", "metadata": {}, "source": [ "#### GTC\n", "\n", "Here we compute the geometric terrain correction.\n", "\n", "Input parameters:\n", "- `product_urlpath`: product path\n", "- `measurement_group`: band to be processed in the form {swath}/{polarization} (see [xarray-sentinel](https://pypi.org/project/xarray-sentinel/) for more details)\n", "- `dem_urlpath`: path to the input DEM. sarsen supports all DEMs supported by GDAL/Proj for ECEF-translation. \n", "- `interp_method`: interpolation method, sarsen supports all interpolation methods supported by [xarray.Dataset.interp](https://xarray.pydata.org/en/stable/generated/xarray.Dataset.interp.html)\n", "- `chunks`: dask chunks\n", "- `output_urlpath`: output path\n", "\n", "The output is the input SAR image resampled on DEM coordinates" ] }, { "cell_type": "code", "execution_count": 15, "id": "41add520", "metadata": {}, "outputs": [], "source": [ "gtc_path = os.path.basename(product_folder) + \".GTC.tif\"\n", "\n", "gtc = apps.terrain_correction(\n", " grd_local_path,\n", " measurement_group,\n", " dem_path,\n", " interp_method=\"nearest\",\n", " override_product_files=\"{dirname}/{prefix}{swath}-{polarization}{ext}\",\n", " chunks=2048,\n", " output_urlpath=gtc_path,\n", ")" ] }, { "cell_type": "code", "execution_count": 16, "id": "27a84435", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 447, "width": 629 } }, "output_type": "display_data" } ], "source": [ "_ = gtc.plot(vmax=0.4)" ] }, { "cell_type": "markdown", "id": "c7399622", "metadata": {}, "source": [ "#### RTC\n", "`sarsen` implements the radiometric terrain-correction [Gamma Flattening](https://ieeexplore.ieee.org/document/5752845) algorithm.\n", "\n", "#### Input parameters\n", "- `correct_radiometry`: `correct_radiometry`: default `None`. If `correct_radiometry=None`the radiometric terrain correction is not applied. `correct_radiometry=gamma_bilinear` applies the gamma flattening classic algorithm using bilinear interpolation to compute the weights. `correct_radiometry=gamma_nearest` applies the gamma flattening using nearest neighbours instead of bilinear interpolation. 'gamma_nearest' significantly reduces the processing time.\n", "- `grouping_area_factor`: scaling factor for the size of the image pixel where the areas are summed. By default, the `grouping_area_factor` is `(1, 1)`, which corresponds to Sentinel-1 input product pixel size. The `grouping_area_factor` shall be increased if the DEM resolution is lower than the Sentinel-1 input product resolution to avoid gaps and distortions the normalization factor. It can be also used to to speed up the computation or the DEM resolution is lower than the Sentinel-1 input product resolution. \n", "\n", "\n", "**Note**: The `grouping_area_factor` can be increased (i) to speed up the processing or (ii) when the input DEM resolution is low. The Gamma Flattening usually works properly if the pixel size of the input DEM is much smaller than the pixel size of the input Sentinel-1 product. Otherwise, the output may have radiometric distortions. This problem can be avoided by increasing the `grouping_area_factor`. Be aware that `grouping_area_factor` too high may degrade the final result." ] }, { "cell_type": "code", "execution_count": 17, "id": "8416dd32-afa0-434a-8c9c-5d44bdd06cf9", "metadata": {}, "outputs": [], "source": [ "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\", \"divide by zero\", RuntimeWarning)" ] }, { "cell_type": "code", "execution_count": 18, "id": "bfdd3e96", "metadata": {}, "outputs": [], "source": [ "rtc_path = os.path.basename(grd_local_path) + \".RTC.tif\"\n", "\n", "rtc = apps.terrain_correction(\n", " grd_local_path,\n", " measurement_group,\n", " dem_path,\n", " interp_method=\"nearest\",\n", " override_product_files=\"{dirname}/{prefix}{swath}-{polarization}{ext}\",\n", " correct_radiometry=\"gamma_nearest\",\n", " output_urlpath=rtc_path,\n", " grouping_area_factor=(1, 1),\n", ")" ] }, { "cell_type": "code", "execution_count": 19, "id": "849a9019", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 447, "width": 629 } }, "output_type": "display_data" } ], "source": [ "_ = rtc.plot(vmax=0.4)" ] }, { "cell_type": "markdown", "id": "646acd93", "metadata": {}, "source": [ "### Gamma vs Local incidence Angle\n", "\n", "The backscattering depends on the ground material, the roughness and the local incidence angle. Small values of the local incidence angle correspond usually to higher values of backscattering, while high values of radar incidence angle corresponds to zones with low backscattering. The aim of the radiometric terrain correction is to remove the contribution due to the local incidence angle. \n", "For evaluating the goodness of the radiometric correction we will compare the incidence angle versus backscattering before and after the correction." ] }, { "cell_type": "markdown", "id": "3bc92d66", "metadata": {}, "source": [ "#### Compute local incidence angle and ellipsoid incidence angle" ] }, { "cell_type": "markdown", "id": "f93af76d", "metadata": {}, "source": [ "![SNOWFALL](https://forum.step.esa.int/uploads/default/original/3X/e/3/e372bf0efa0133c63977111b15614bc0c7d96fd8.png)" ] }, { "cell_type": "code", "execution_count": 20, "id": "e8499527-c404-4c45-b266-c0ec077db4d9", "metadata": {}, "outputs": [], "source": [ "dem = scene.open_dem_raster(dem_path)\n", "\n", "\n", "orbit_ecef = xr.open_dataset(\n", " grd_local_path,\n", " engine=\"sentinel-1\",\n", " group=f\"{measurement_group}/orbit\",\n", " override_product_files=\"{dirname}/{prefix}{swath}-{polarization}{ext}\",\n", ")\n", "dem_ecef = scene.convert_to_dem_ecef(dem, source_crs=t_srs)\n", "dem_ecef_1 = scene.convert_to_dem_ecef(dem + 1, source_crs=t_srs)\n", "ellissoid_normal = dem_ecef - dem_ecef_1" ] }, { "cell_type": "code", "execution_count": 21, "id": "a5f3dd32-007b-4c1c-9cda-88592fe78dd4", "metadata": {}, "outputs": [], "source": [ "acquisition = apps.simulate_acquisition(\n", " dem_ecef, orbit_ecef.position, coordinate_conversion=None, correct_radiometry=True\n", ")\n", "oriented_area = scene.compute_dem_oriented_area(dem_ecef)\n", "dem_normal = -oriented_area / np.sqrt(xr.dot(oriented_area, oriented_area, dims=\"axis\"))\n", "\n", "# compute direction target sensor\n", "orbit_interpolator = orbit.OrbitPolyfitIterpolator.from_position(orbit_ecef.position)\n", "position_ecef = orbit_interpolator.position()\n", "velocity_ecef = orbit_interpolator.velocity()\n", "acquisition = geocoding.backward_geocode(dem_ecef, orbit_ecef.position, velocity_ecef)\n", "slant_range = np.sqrt((acquisition.dem_distance**2).sum(dim=\"axis\"))\n", "dem_direction = acquisition.dem_distance / slant_range" ] }, { "cell_type": "code", "execution_count": 22, "id": "920ee292-f695-44d4-9178-950fb5703a4b", "metadata": {}, "outputs": [], "source": [ "angle = np.arccos(xr.dot(dem_normal, dem_direction, dims=\"axis\"))" ] }, { "cell_type": "code", "execution_count": 23, "id": "a2797f82-5471-4a75-9878-317e9b23865d", "metadata": {}, "outputs": [], "source": [ "ellipsoid_incidence = np.arccos(xr.dot(ellissoid_normal, dem_direction, dims=\"axis\"))" ] }, { "cell_type": "markdown", "id": "5c69085f", "metadata": {}, "source": [ "#### Load Data" ] }, { "cell_type": "code", "execution_count": 24, "id": "a45cd942", "metadata": {}, "outputs": [], "source": [ "ellipsoid_incidence = ellipsoid_incidence.load()\n", "angle = angle.load()\n", "gtc = gtc.load()\n", "rtc = rtc.load()" ] }, { "cell_type": "markdown", "id": "5311f17d-aa35-4725-8e28-dd91719cf91c", "metadata": {}, "source": [ "#### GTC and RTC distribution of foreslope, backslope, and flat terrain \n", "\n", "We classify taking into account the local incidence angle and the ellipsoid incidence angle\n", " the image in three zone:\n", "- foreslope, small incidence angles that correspond to high backscattering zones.\n", "- flat terrain\n", "- backslope, big incidence angles that correspond to low backscattering zones\n", "\n", "\n", "We will show that the radiometric correction mitigate the effect of the foreslope and backslope." ] }, { "cell_type": "code", "execution_count": 25, "id": "fec9e2e8-0dc7-49e5-a9ea-a12bbce248e9", "metadata": {}, "outputs": [], "source": [ "# foreslope backslope and flat terrain\n", "gtc_foreslope = xr.where(angle < ellipsoid_incidence - 0.20, gtc, np.nan)\n", "gtc_backslope = xr.where(angle > ellipsoid_incidence + 0.20, gtc, np.nan)\n", "gtc_flat = xr.where(\n", " (angle < ellipsoid_incidence + 0.20) & (angle > ellipsoid_incidence - 0.20),\n", " gtc,\n", " np.nan,\n", ")\n", "\n", "rtc_foreslope = xr.where(angle < ellipsoid_incidence - 0.20, rtc, np.nan)\n", "rtc_backslope = xr.where(angle > ellipsoid_incidence + 0.20, rtc, np.nan)\n", "rtc_flat = xr.where(\n", " (angle < ellipsoid_incidence + 0.20) & (angle > ellipsoid_incidence - 0.20),\n", " rtc,\n", " np.nan,\n", ")" ] }, { "cell_type": "code", "execution_count": 26, "id": "e3abd619-34e2-4fec-845d-89e85adf4049", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 631, "width": 2109 } }, "output_type": "display_data" } ], "source": [ "f, axes = plt.subplots(nrows=1, ncols=3, figsize=(30, 9))\n", "\n", "dem_raster.plot(ax=axes[0])\n", "axes[0].grid(c=\"red\")\n", "axes[0].set_title(\"DEM\")\n", "\n", "gtc.plot(ax=axes[1], vmax=0.4)\n", "axes[1].grid(c=\"red\")\n", "axes[1].set_title(\"GTC\")\n", "\n", "\n", "angle.plot(ax=axes[2])\n", "axes[2].grid(c=\"red\")\n", "axes[2].set_title(\"Incidence Angle\")\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": 27, "id": "8923dfd6-b740-4a72-9f88-a67e07d47df6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 631, "width": 2114 } }, "output_type": "display_data" } ], "source": [ "f, axes = plt.subplots(nrows=1, ncols=3, figsize=(30, 9))\n", "\n", "gtc_backslope.plot(ax=axes[0], vmax=0.4)\n", "axes[0].grid(c=\"red\")\n", "axes[0].set_title(\"GTC basckslope\")\n", "\n", "gtc_foreslope.plot(ax=axes[1], vmax=0.4)\n", "axes[1].grid(c=\"red\")\n", "axes[1].set_title(\"GTC foreslope\")\n", "\n", "\n", "gtc_flat.plot(ax=axes[2], vmax=0.4)\n", "axes[2].grid(c=\"red\")\n", "axes[2].set_title(\"GTC flat terrain\")\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "id": "10bcc74a-902c-4fb3-ac85-f0012e307f68", "metadata": {}, "source": [ "#### GTC and RTC distribution of foreslope backslope and flat terrain in db" ] }, { "cell_type": "code", "execution_count": 28, "id": "12ef90a7", "metadata": {}, "outputs": [], "source": [ "gtc_backslope_db = np.log(gtc_backslope)\n", "gtc_foreslope_db = np.log(gtc_foreslope)\n", "gtc_flat_db = np.log(gtc_flat)\n", "\n", "rtc_backslope_db = np.log(rtc_backslope)\n", "rtc_foreslope_db = np.log(rtc_foreslope)\n", "rtc_flat_db = np.log(rtc_flat)" ] }, { "cell_type": "code", "execution_count": 29, "id": "68081745-5198-441f-87bf-b9b7e49d81ed", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 0, 'gamma_T DB')" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 724, "width": 1717 } }, "output_type": "display_data" } ], "source": [ "f, axes = plt.subplots(nrows=1, ncols=2, figsize=(30, 12))\n", "\n", "gtc_backslope_db_h = xr.plot.hist(\n", " gtc_backslope_db,\n", " range=[-10, 5],\n", " bins=50,\n", " color=[0, 1, 0, 0.4],\n", " density=True,\n", " ax=axes[0],\n", " label=\"backslope\",\n", ")\n", "gtc_foreslope_db_h = xr.plot.hist(\n", " gtc_foreslope_db,\n", " range=[-10, 5],\n", " bins=50,\n", " color=[1, 0, 0, 0.4],\n", " density=True,\n", " ax=axes[0],\n", " label=\"foreslope\",\n", ")\n", "gtc_flat_db_h = xr.plot.hist(\n", " gtc_flat_db,\n", " range=[-10, 5],\n", " bins=50,\n", " color=[0, 0, 1, 0.4],\n", " density=True,\n", " ax=axes[0],\n", " label=\"flat\",\n", ")\n", "\n", "rtc_backslope_db_h = xr.plot.hist(\n", " rtc_backslope_db,\n", " range=[-10, 5],\n", " bins=50,\n", " color=[0, 1, 0, 0.4],\n", " density=True,\n", " ax=axes[1],\n", " label=\"backslope\",\n", ")\n", "rtc_foreslope_db_h = xr.plot.hist(\n", " rtc_foreslope_db,\n", " range=[-10, 5],\n", " bins=50,\n", " color=[1, 0, 0, 0.4],\n", " density=True,\n", " ax=axes[1],\n", " label=\"foreslope\",\n", ")\n", "rtc_flat_db_h = xr.plot.hist(\n", " rtc_flat_db,\n", " range=[-10, 5],\n", " bins=50,\n", " color=[0, 0, 1, 0.4],\n", " density=True,\n", " ax=axes[1],\n", " label=\"flat\",\n", ")\n", "axes[0].legend()\n", "axes[1].legend()\n", "\n", "axes[0].set_title(\"GTC distribution\")\n", "axes[0].set_xlabel(\"beta DB\")\n", "\n", "axes[1].set_title(\"RTC distribution\")\n", "axes[1].set_xlabel(\"gamma_T DB\")" ] }, { "cell_type": "markdown", "id": "cff38443-b393-4ace-bf3f-bc64318347a1", "metadata": {}, "source": [ "**Before the radiometic correction the distribution of backslope flat terrain and foreslope are separate, after they overlap completely.**" ] }, { "cell_type": "markdown", "id": "456b4611", "metadata": {}, "source": [ "#### RTC and GTC scatterplots incidence angle vs backcattering\n", "\n", "We compare the incidence_angle-backscattering scatterplots of the GTC and the RTC to show that in the former there is a correlation between the incidence angle and backscattering, while in the RTC the correlation has been removed." ] }, { "cell_type": "code", "execution_count": 30, "id": "82786f69", "metadata": {}, "outputs": [], "source": [ "def makeColours(vals):\n", " colours = np.zeros((len(vals), 3))\n", " norm = Normalize(vmin=vals.min(), vmax=vals.max())\n", "\n", " # Can put any colormap you like here.\n", " colours = [cm.ScalarMappable(norm=norm, cmap=\"plasma\").to_rgba(val) for val in vals]\n", "\n", " return colours\n", "\n", "\n", "def plot_density(x, y, ax=None):\n", " if ax is None:\n", " ax = plt\n", " mask = np.logical_and(np.isfinite(y), np.isfinite(x))\n", " samples = np.stack([x[mask], y[mask]])\n", "\n", " densObj = kde(samples)\n", " colours = makeColours(densObj.evaluate(samples))\n", " fig = ax.scatter(samples[0], samples[1], color=colours)\n", " return fig" ] }, { "cell_type": "code", "execution_count": 31, "id": "92dc79c9-099b-4c28-b41d-a46ce68cdb87", "metadata": {}, "outputs": [], "source": [ "gtc_db = np.log(gtc).load()\n", "rtc_db = np.log(rtc).load()\n", "angle = angle.load()" ] }, { "cell_type": "code", "execution_count": 32, "id": "498a010d-6465-4cce-89c9-1feeb5cd5a0d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 0, 'local incidence angle')" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 724, "width": 1734 } }, "output_type": "display_data" } ], "source": [ "f, axes = plt.subplots(nrows=1, ncols=2, figsize=(30, 12))\n", "\n", "fig = plot_density(\n", " angle[::10, ::10].data.ravel(), gtc_db[::10, ::10].data.ravel(), ax=axes[0]\n", ")\n", "fig = plot_density(\n", " angle[::10, ::10].data.ravel(), rtc_db[::10, ::10].data.ravel(), ax=axes[1]\n", ")\n", "\n", "axes[0].set_title(\"GTC\")\n", "axes[0].set_ylabel(\"beta DB\")\n", "axes[0].set_xlabel(\"local incidence angle\")\n", "\n", "axes[1].set_title(\"RTC\")\n", "axes[1].set_ylabel(\"gamma_T DB\")\n", "axes[1].set_xlabel(\"local incidence angle\")" ] }, { "cell_type": "markdown", "id": "e9a66c4d-ec35-49a9-8b71-1594615d6403", "metadata": {}, "source": [ "**Note that in the first plot there is a clear correlation between the angle of incidence and backscattering, while in the second graph there is no correlation.**" ] } ], "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": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }