Surface Water Dynamics - Time Series Analysis with Sentinel-1

By: @WalidGharianiEAGLE | DHI

Introduction

Water is a vital part of Earth ecosystems and life, supporting biodiversity, and sustaining human livelihoods. Monitoring surface water dynamics (occurence, frequency and change) is important for managing resources, and mitigating natural hazards such as floods and droughts. Wetlands in particular are unique ecosystems under the influence of precipitaion, hydrological processes and coastal dynamics, which contribute in shaping the habitats, and species diversity. Within this context, the Keta and Songor Ramsar Sites in southeastern Ghana present a dynamic coastal wetland system where water levels fluctuate due to rainfall, tidal influence, and lagoon-river interactions.

These ramsar sites consists of different wetlands classes such as marshes, floodplains, mangroves, and seasonally inundated grasslands. These unique sites are also critical habitats for thousands of resident and migratory birds, fish, and sea turtles, while supporting local livelihoods through fishing, agriculture, and salt production. Threfore an effective monitoring of water dynamics of these ecosystems is essential for conserving biodiversity and sustaining community resources.

Sentinel-1 Synthetic Aperture Radar (SAR) can detect water under all weather conditions, overcoming limitations of optical imagery caused by cloud cover. By processing and analyzing Sentinel-1 time series and its radar backscatter (VH, VV), water occurrence, frequency and dynamics in Ramsar wetlands could be detected and monitored providing valuable insights of these complex wetlands sites.

In this notebook, we will explore how to monitor surface water dynamics in coastal wetlands using Sentinel-1 time series. We will process radar backscatter data to detect water occurrence, frequency, and change over time.

What we will learn

  • 🛰️ Accessing Sentinel-1 Ground Range Detected (GRD) .zarr data via the EODC STAC API using pystac-client.
  • 🛠️ Preprocessing radar imagery including spatial subsetting, radiometric calibration, speckle filtering, georeferencing, and regridding.
  • 🌊 Generating surface water masks using an adaptative thresholding algorithm.
  • 📊 Time series analysis of water dynamics.

Import libraries

import datetime as dt
from tqdm import tqdm

import fsspec
import numpy as np
import xarray as xr
import hvplot.xarray
import pandas as pd
import matplotlib.pyplot as plt
from pystac_client import Client

import session_info
import warnings

warnings.filterwarnings("ignore", category=UserWarning, module="pkg_resources")

Helper functions

Before starting the analysis workflow, we import several helper functions from zarr_s1_utils.py that implement key processing steps for Sentinel-1.

from zarr_s1_utils import (
    subset,
    radiometric_calibration,
    lee_filter_dask,
    regrid,
    xr_threshold_otsu,
)
session_info.show()
Click to view session information 
-----
fsspec              2026.2.0
hvplot              0.12.2
matplotlib          3.10.8
numpy               2.4.4
pandas              2.3.3
pystac_client       0.8.6
session_info        v1.0.1
tqdm                4.67.3
xarray              2026.2.0
zarr_s1_utils       NA
-----
Click to view modules imported as dependencies 
81d243bd2c585b0f4821__mypyc NA
PIL                         12.2.0
affine                      2.4.0
annotated_types             0.7.0
anyio                       NA
arrow                       1.4.0
asttokens                   NA
attr                        26.1.0
attrs                       26.1.0
babel                       2.18.0
bleach                      6.3.0
bokeh                       3.9.0
boto3                       1.42.84
botocore                    1.42.84
bs4                         4.14.3
cachetools                  7.0.5
certifi                     2026.02.25
charset_normalizer          3.4.7
click                       8.3.2
cloudpickle                 3.1.2
colorama                    0.4.6
colorcet                    3.1.0
comm                        0.2.3
cycler                      0.12.1
cython_runtime              NA
dask                        2026.3.0
dateutil                    2.9.0.post0
debugpy                     1.8.20
decorator                   5.2.1
defusedxml                  0.7.1
dotenv                      NA
executing                   2.2.1
fastjsonschema              NA
fqdn                        NA
func_timeout                4.3.5
geopandas                   1.1.3
h3                          4.4.2
holoviews                   1.22.1
idna                        3.11
igraph                      0.11.9
ipykernel                   7.2.0
isoduration                 NA
jedi                        0.19.2
jinja2                      3.1.6
jmespath                    1.1.0
joblib                      1.5.3
json5                       0.14.0
jsonpointer                 3.1.1
jsonschema                  4.26.0
jsonschema_specifications   NA
jupyter_events              0.12.0
jupyter_server              2.17.0
jupyterlab_pygments         0.3.0
jupyterlab_server           2.28.0
kiwisolver                  1.5.0
lark                        1.3.1
lazy_loader                 0.5
markupsafe                  3.0.3
matplotlib_inline           0.2.1
mistune                     3.2.0
mpl_toolkits                NA
msgpack                     1.1.2
narwhals                    2.19.0
nbclient                    0.10.4
nbconvert                   7.17.1
nbformat                    5.10.4
networkx                    3.6.1
numexpr                     2.14.1
odc                         NA
osmnx                       2.0.6
packaging                   26.0
pandocfilters               NA
panel                       1.8.10
param                       2.3.3
parso                       0.8.6
planetary_computer          1.0.0
platformdirs                4.9.6
prometheus_client           NA
prompt_toolkit              3.0.52
psutil                      7.2.2
pure_eval                   0.2.3
pyarrow                     23.0.1
pydantic                    2.12.5
pydantic_core               2.41.5
pydev_ipython               NA
pydevconsole                NA
pydevd                      3.2.3
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.20.0
pymannkendall               1.4.3
pyparsing                   3.3.2
pyproj                      3.7.2
pystac                      1.14.3
pythonjsonlogger            NA
pytz                        2026.1.post1
pyviz_comms                 3.0.6
rasterio                    1.5.0
referencing                 NA
regex                       2026.4.4
requests                    2.33.1
rfc3339_validator           0.1.4
rfc3986_validator           0.1.1
rfc3987_syntax              NA
rioxarray                   0.22.0
rpds                        NA
s3transfer                  0.16.0
scipy                       1.17.1
send2trash                  NA
shapely                     2.1.2
six                         1.17.0
skimage                     0.26.0
sklearn                     1.8.0
soupsieve                   2.8.3
stack_data                  0.6.3
tblib                       3.2.2
texttable                   1.7.0
threadpoolctl               3.6.0
tinycss2                    1.4.0
tlz                         1.1.0
toolz                       1.1.0
tornado                     6.5.5
traitlets                   5.14.3
typing_extensions           NA
typing_inspection           NA
uri_template                NA
urllib3                     2.6.3
wcwidth                     0.6.0
webcolors                   NA
webencodings                0.5.1
websocket                   1.9.0
xxhash                      NA
yaml                        6.0.3
zmq                         27.1.0
zoneinfo                    NA

-----
IPython             9.12.0
jupyter_client      8.8.0
jupyter_core        5.9.1
jupyterlab          4.5.6
notebook            7.5.5
-----
Python 3.12.13 | packaged by conda-forge | (main, Mar  5 2026, 16:50:00) [GCC 14.3.0]
Linux-6.17.0-1010-azure-x86_64-with-glibc2.39
-----
Session information updated at 2026-04-16 09:05

Preprocessing

Spatial subset

To focus our analysis over the chosen AOI, we can efficiently crop our dataset using a spatial subset. The function subset() determine the slices in azimuth_time and ground_range that cover the AOI, and then extract and mask the corresponding portion of the GRD dataset.

The function subset() takes the following keyword arguments:

  • grd:The GRD dataset to be cropped (the radar image in azimuth and ground range coordinates).
  • gcp_ds: The GCP (Ground Control Points) dataset containing the latitude and longitude grids used to geolocate the GRD image.
  • aoi_bounds: The geographic bounding box of the Area of Interest, given as: [min_lon, min_lat, max_lon, max_lat].
  • offset: The number of GCP grid cells to include around the AOI center. This adds a small margin around the AOI to ensure the cropped region fully covers it.
grd_subset = subset(grd=grd, gcp_ds=gcp, aoi_bounds=aoi_bounds, offset=1)
grd_subset
<xarray.DataArray 'grd' (azimuth_time: 1061, ground_range: 1316)> Size: 6MB
dask.array<where, shape=(1061, 1316), dtype=float32, chunksize=(1061, 1316), chunktype=numpy.ndarray>
Coordinates:
  * azimuth_time  (azimuth_time) datetime64[ns] 8kB 2024-12-28T18:10:22.94533...
    line          (azimuth_time) float64 8kB 3.356e+03 3.357e+03 ... 4.416e+03
  * ground_range  (ground_range) float64 11kB 3.683e+04 3.684e+04 ... 4.998e+04
    pixel         (ground_range) float64 11kB 3.683e+03 3.684e+03 ... 4.998e+03
    latitude      (azimuth_time, ground_range) float64 11MB 5.739 ... 5.859
    longitude     (azimuth_time, ground_range) float64 11MB 0.6134 ... 0.7102
Attributes:
    _eopf_attrs:  {'coordinates': ['azimuth_time', 'line', 'pixel', 'ground_r...
    dtype:        <u2
    long_name:    measurement data set for GRD IW
grd_subset.plot(robust=True, cmap="cividis")
plt.show()

Radiometric calibration

In order to get a meanigful physical properties of features in the SAR scene that could be used for quantitative analysis, we need to apply a radiometric calibration on the backscatter values. This step converts the backscatter into a calibrated normalized radar cross section, correcting for incidence angle and sensor characteristics ensuring SAR images from different dates or viewing geometries are directly comparable.

  • Reference: https://step.esa.int/docs/tutorials/S1TBX%20SAR%20Basics%20Tutorial.pdf

The radiometric calibration is done using the radiometric_calibration() function which takes the following keyword arguments:

  • grd:The GRD data array to be calibrated.
  • calibration_ds: The calibration dataset containing the radiometric calibration lookup tables provided with the product. This dataset is interpolated to the GRD grid before being applied.
  • calibration_type: The name of the calibration parameter to use from the calibration dataset. In this case, sigma_nought is used to compute the sigma nought backscatter coefficient.

To have a look into the available data vars within the calibration dataset that could be used for the radiometric calibration

calibration.data_vars
Data variables:
    beta_nought   (azimuth_time, ground_range) float32 68kB dask.array<chunksize=(27, 632), meta=np.ndarray>
    dn            (azimuth_time, ground_range) float32 68kB dask.array<chunksize=(27, 632), meta=np.ndarray>
    gamma         (azimuth_time, ground_range) float32 68kB dask.array<chunksize=(27, 632), meta=np.ndarray>
    sigma_nought  (azimuth_time, ground_range) float32 68kB dask.array<chunksize=(27, 632), meta=np.ndarray>

We apply the calibration, obtaining:

sigma_0 = radiometric_calibration(
    grd=grd_subset, calibration_ds=calibration, calibration_type="sigma_nought"
)
sigma_0
<xarray.DataArray (azimuth_time: 1061, ground_range: 1316)> Size: 6MB
dask.array<pow, shape=(1061, 1316), dtype=float32, chunksize=(1061, 1316), chunktype=numpy.ndarray>
Coordinates:
  * azimuth_time  (azimuth_time) datetime64[ns] 8kB 2024-12-28T18:10:22.94533...
  * ground_range  (ground_range) float64 11kB 3.683e+04 3.684e+04 ... 4.998e+04
    latitude      (azimuth_time, ground_range) float64 11MB 5.739 ... 5.859
    longitude     (azimuth_time, ground_range) float64 11MB 0.6134 ... 0.7102
sigma_0.plot(robust=True, cmap="cividis", cbar_kwargs={"label": "Sigma Nought"})

Speckel filtering

Raw SAR imagery is characterized by “grainy” or “salt and pepper” effect caused by random constructive and destructive interference, known as speckle. In order to reduce this effect and noise, we apply the spatial Lee Filter (Lee et al., 2009) that averages the pixel values while preserving edges.

For speckel filtering we will use the lee_filter_dask() function which takes the following keyword arguments:

  • da: The input xarray.DataArray to be filtered. Here, sigma_0 is the radiometrically calibrated GRD subset.
  • size: The size of the square moving window used to compute the local statistics for the Lee filter. Odd numbers are recommended (e.g., 5×5). Larger windows produce stronger smoothing but may reduce spatial detail.
sigma_0_spk = lee_filter_dask(da=sigma_0, size=5)
sigma_0_spk
<xarray.DataArray (azimuth_time: 1061, ground_range: 1316)> Size: 6MB
dask.array<transpose, shape=(1061, 1316), dtype=float32, chunksize=(1061, 1316), chunktype=numpy.ndarray>
Coordinates:
  * azimuth_time  (azimuth_time) datetime64[ns] 8kB 2024-12-28T18:10:22.94533...
  * ground_range  (ground_range) float64 11kB 3.683e+04 3.684e+04 ... 4.998e+04
    latitude      (azimuth_time, ground_range) float64 11MB 5.739 ... 5.859
    longitude     (azimuth_time, ground_range) float64 11MB 0.6134 ... 0.7102
Attributes:
    attrs:    {'speckle filter method': 'lee'}

Obtaining a cleaner image

sigma_0_spk.plot(robust=True, cmap="cividis", cbar_kwargs={"label": "Sigma Nought"})
plt.title("Sigma Nought after Lee Filter")
plt.show()

Georefrecing and regredding

Sentinel-1 GRD Zarr comes in irregular image geometry that does not align with common geographic grids. In order to conduct spatial analysis, and scenes comparison, and mapping, we need to georefrence and resample the data onto a regular grid.
We will use an ODC GeoBox along with with scipy.griddata to interpolate the SAR values onto a consistent latitude-longitude grid. These utilities are available within the regrid() function which takes the following keyword arguments:

  • da: The input xarray.DataArray to regrid. Here, sigma_0_spk is the speckle-filtered GRD subset.
  • bounds: The geographic bounding box for the output grid: [min_lon, min_lat, max_lon, max_lat].
  • resolution: Tuple (dx, dy) defining the grid spacing in the coordinate units of the CRS. For example, 10 / 111320 degrees corresponds roughly to 10 meters.
  • crs: Coordinate reference system of the output grid (default is "EPSG:4326" for WGS84).
  • method: Interpolation method for mapping irregularly spaced data to the regular grid. Options include "nearest", "linear", or "cubic".
resolution = 10 / 111320  # approx 10 meters in degrees
crs = "epsg:4326"

sigma_0_spk_geo = regrid(
    da=sigma_0_spk,
    bounds=aoi_bounds,
    resolution=resolution,
    crs=crs,
    method="nearest",  # "linear" or "cubic"
)
sigma_0_spk_geo
<xarray.DataArray (y: 868, x: 1175)> Size: 8MB
array([[0.02022304, 0.02022304, 0.0271727 , ..., 0.00480022, 0.00438391,
        0.00438391],
       [0.02022304, 0.02022304, 0.0271727 , ..., 0.00480022, 0.0070433 ,
        0.0070433 ],
       [0.02168461, 0.02168461, 0.02967146, ..., 0.00716109, 0.00948818,
        0.00948818],
       ...,
       [0.00227514, 0.00227514, 0.00214053, ..., 0.00365857, 0.00334549,
        0.00260707],
       [0.00176132, 0.00176132, 0.00177638, ..., 0.00396993, 0.00376833,
        0.00376833],
       [0.00176132, 0.00177638, 0.00177638, ..., 0.00396993, 0.00376833,
        0.00376833]], shape=(868, 1175))
Coordinates:
  * y            (y) float64 7kB 5.838 5.838 5.838 5.837 ... 5.76 5.76 5.76 5.76
  * x            (x) float64 9kB 0.6091 0.6092 0.6093 ... 0.7144 0.7145 0.7146
    spatial_ref  int64 8B 0
Attributes:
    interpolation:  nearest
    odcgeobox:      GeoBox((868, 1175), Affine(8.983111749910169e-05, 0.0, 0....
sigma_0_spk_geo.odcgeobox

GeoBox

Dimensions
1,175x868
EPSG
4326
Resolution
8.98311e-05°
Cell
200px
WKT
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
sigma_0_spk_geo.plot(robust=True, cmap="cividis", cbar_kwargs={"label": "Sigma Nought"})
plt.title("Regridded Sigma Nought after Lee Filter")
plt.show()

Convert backscatter to dB

We convert the regridded sigma_0 backscatter intensity from Linear scale to decibels (dB) using a logarithmic transformation. This enhances contrast and simplifies statistical analysis and interpretation of the image. It is considered a standard approach for representing SAR intensity.

sigma_0_spk_geo_db = 10 * np.log10(sigma_0_spk_geo)
sigma_0_spk_geo_db
<xarray.DataArray (y: 868, x: 1175)> Size: 8MB
array([[-16.94153556, -16.94153556, -15.65867213, ..., -23.18738561,
        -23.58138141, -23.58138141],
       [-16.94153556, -16.94153556, -15.65867213, ..., -23.18738561,
        -21.52223657, -21.52223657],
       [-16.63848472, -16.63848472, -15.27661105, ..., -21.45021141,
        -20.22816859, -20.22816859],
       ...,
       [-26.42991425, -26.42991425, -26.69479123, ..., -24.36688424,
        -24.75540359, -25.8384678 ],
       [-27.54160746, -27.54160746, -27.50464296, ..., -24.01217139,
        -24.2385137 , -24.2385137 ],
       [-27.54160746, -27.50464296, -27.50464296, ..., -24.01217139,
        -24.2385137 , -24.2385137 ]], shape=(868, 1175))
Coordinates:
  * y            (y) float64 7kB 5.838 5.838 5.838 5.837 ... 5.76 5.76 5.76 5.76
  * x            (x) float64 9kB 0.6091 0.6092 0.6093 ... 0.7144 0.7145 0.7146
    spatial_ref  int64 8B 0
Attributes:
    interpolation:  nearest
    odcgeobox:      GeoBox((868, 1175), Affine(8.983111749910169e-05, 0.0, 0....
sigma_0_spk_geo_db.plot(
    robust=True, cmap="cividis", cbar_kwargs={"label": "Sigma Nought dB"}
)
plt.title("Calibrated Sigma Nought in dB")
plt.show()

sigma_0_spk_geo_db.hvplot.image(
    x="x",
    y="y",
    robust=True,
    cmap="cividis",
    title="SAR GRD",
)