Using Sentinel-2 images to monitor mouth of the river Syr Darya and Aral Lake and analyzing it with NDWI index on Creodias

The Syr Darya supplies the North Aral Sea, where the availability of water matters to aquatic habitats and fisheries. Near the river mouth, a shift in the visible shoreline or the loss of a connection between water bodies can point to a change worth investigating. Satellite observations provide a way to compare these features over time and identify places that need closer inspection. This example uses three Sentinel-2 acquisitions to produce comparable NDWI maps, helping you look for possible changes in surface-water coverage and decide where additional imagery, river-flow records or field observations are needed.

Sentinel-2A Level-2A images from Creodias show water-related spectral patterns near the mouth of the Syr Darya river and Aral Lake in Kazakhstan, from May through July 2023. The example downloads three complete satellite products to the local disk, creates Normalized Difference Water Index (NDWI) rasters from green and near-infrared surface reflectance, and compares acquisitions from 26 May, 5 June and 5 July 2023. NDWI does not measure water depth or volume, and the available acquisitions do not provide a daily record.

What is NDWI?

NDWI is a normalized difference index calculated based on the reflectance of light in two different spectral bands: green (G) and near-infrared (NIR). The formula for calculating NDWI is as follows:

where:

• G - reflectance in the green band,

• NIR - reflectance in the near-infrared band.

NDWI values range from -1 to 1, where higher values indicate the presence of water, while the lower values may indicate other land cover types, such as vegetation or urban areas.

Advantages of using NDWI indicator

Raster images containing the NDWI indicator can be used in various fields, including:

Water Resource Monitoring

Comparing water extent in lakes, rivers and reservoirs, including areas affected by drought or excessive rainfall.

Crisis Management

Rapid identification of flooded areas during floods and assessment of damages.

Environmental Protection

Monitoring wetlands and other aquatic ecosystems to preserve their biodiversity.

Urban Planning

Analyzing the impact of urban development on nearby water resources.

Agriculture

Assessing irrigation of agricultural fields and monitoring plant health.

Prerequisites

No. 1 Account and product downloads

A Creodias account and an access token are required for the product downloads. The public catalogue search returns metadata only. The authentication procedure in Using curl and wget to download EODATA products from Creodias leaves the token in the Bash variable ACCESS_TOKEN. The download command runs in that same terminal, which needs Python and curl. Credentials and tokens are not entered into notebook cells or saved in the notebook.

No. 2 Local Jupyter environment and disk space

The example runs in a local Jupyter installation. The Jupyter installation instructions describe its setup. Satellite products are downloaded and extracted into Downloads/aral-data inside the current user’s home directory. No EODATA filesystem mount is needed. The three archives total about 3.4 GB; at least 9 GiB of free space is required for the archives, extracted products and output rasters.

The files for this example are sentinel2_aral_lake_branded.ipynb and download-aral-comparison.py. Both belong in the local Downloads folder. The notebook has the Creodias endpoint already configured.

The kernel needs requests, NumPy, Rasterio and Matplotlib. Rasterio must support JPEG2000 to read the Sentinel-2 bands. The following command runs in a Jupyter Code cell with Shift+Enter; it installs packages into the environment used by that kernel.

%pip install requests numpy rasterio matplotlib

Installation makes the packages available; the import cells load them into the current Python session. After a kernel restart, the parameter, import and function-definition cells must run again before processing. Markdown cells contain explanations and are not Python commands.

No. 3 Catalogue and product metadata

The query syntax is described in EOData Catalogue API Manual on Creodias. The complete Level-2A SAFE product contains MTD_MSIL2A.xml, the 10 m B03 and B08 bands, and the 20 m scene classification layer (SCL). The Sentinel-2 processing documentation explains the reflectance offsets and quantification values read from the metadata.

import requests

Define variables to hold the data

Select the acquisition period and area around the Syr Darya river mouth and Aral Lake. The coordinates are defined through BBOX (Bounding Box), as a list of coordinates – Xmin, Ymin, Xmax, Ymax. All coordinates will be defined in EPSG:4326.

In python, we define three variables:

# Timerange of data that we want to recieve
start_date = '2023-05-01'
end_date = '2023-07-31'
# Geometry in form of a BBOX
bbox = [60.2347,46.1731,61.1160,46.6541]

Work with the OData API

The OData API searches EODATA products using spatial and temporal filters. Its response contains each product path in the S3Path field. These remote paths are metadata, not local folders. The download step below supplies the local files used for processing.

Configure the OData endpoint

The Creodias catalogue endpoint is https://datahub.creodias.eu/odata/v1. The notebook and download helper use this endpoint.

ODATA_CATALOGUE_URL = "https://datahub.creodias.eu/odata/v1"
if not ODATA_CATALOGUE_URL.startswith("https://"):
    raise ValueError("Set ODATA_CATALOGUE_URL to the HTTPS OData base URL shown above.")
API_URL = ODATA_CATALOGUE_URL.rstrip("/") + "/Products"

Create the OData query

The query selects Sentinel-2A Level-2A products intersecting the bounding box from 1 May through 31 July 2023, with product cloud cover of at most 10 percent. Processing later crops each granule to the bounding box and applies a pixel-level quality mask. requests encodes the query parameters for the URL.

minx, miny, maxx, maxy = bbox
polygon = (
    f"POLYGON(({minx} {miny},{maxx} {miny},{maxx} {maxy},"
    f"{minx} {maxy},{minx} {miny}))"
)
query_params = {
    "$filter": (
        "Collection/Name eq 'SENTINEL-2' and startswith(Name,'S2A') and "
        "Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' "
        "and att/OData.CSC.StringAttribute/Value eq 'S2MSI2A') and "
        "Attributes/OData.CSC.DoubleAttribute/any(att:att/Name eq 'cloudCover' "
        "and att/OData.CSC.DoubleAttribute/Value le 10) and "
        f"ContentDate/Start ge {start_date}T00:00:00.000Z and "
        f"ContentDate/Start le {end_date}T23:59:59.999Z and "
        f"OData.CSC.Intersects(area=geography'SRID=4326;{polygon}')"
    ),
    "$orderby": "ContentDate/Start asc",
    "$top": 1000,
    "$select": "Id,Name,S3Path",
}
specified_url = requests.Request("GET", API_URL, params=query_params).prepare().url
specified_url

Send the request and follow @odata.nextLink while further result pages are available:

products = []
next_url = specified_url
while next_url:
    response = requests.get(next_url, timeout=90)
    response.raise_for_status()
    page = response.json()
    products.extend(page["value"])
    next_url = page.get("@odata.nextLink")

if not products:
    raise RuntimeError("No Sentinel-2 Level-2A products matched this region and date range. Check catalogue coverage before changing the query parameters.")

Display the first product:

products[0]

Example of the returned fields:

{
  "Id": "49ec07df-edc5-4410-80ee-28312a787128",
  "Name": "S2A_MSIL2A_20230503T064621_N0510_R020_T40TGS_20240903T161656.SAFE",
  "S3Path": "/eodata/Sentinel-2/MSI/L2A_N0500/2023/05/03/S2A_MSIL2A_20230503T064621_N0510_R020_T40TGS_20240903T161656.SAFE"
}

Inspect catalogue paths

The S3Path values describe locations in the remote repository. They do not establish local access to the imagery and are not the processing input for this example.

catalogue_paths = [product["S3Path"] for product in products]
catalogue_paths

Download three acquisitions to the local disk

The comparison in this example requires three separate satellite acquisitions, taken on 26 May, 5 June and 5 July 2023. Each acquisition records the same study area at a different time. The three resulting NDWI maps allow the state of the area to be compared across May, June and July. All three complete products are needed for the three-panel comparison; each comes from tile T40TGS and relative orbit R063.

Acquisitions used for the comparison

Acquisition date	Tile	Relative orbit
26 May 2023	T40TGS	R063
5 June 2023	T40TGS	R063
5 July 2023	T40TGS	R063

All three products are downloaded and extracted into the same local input folder: $HOME/Downloads/aral-data/. Here, $HOME refers to the current user’s home directory. The helper creates aral-data automatically. Its name matters because the notebook later reads comparison-products.json from that exact folder.

A single run of the helper downloads all three products in sequence, extracts their ZIP archives and writes the manifest of local paths. The following command runs in the Bash terminal where ACCESS_TOKEN is set.

printf '%s' "$ACCESS_TOKEN" | python "$HOME/Downloads/download-aral-comparison.py" --catalogue-url "https://datahub.creodias.eu/odata/v1" --data-dir "$HOME/Downloads/aral-data"

After checking the required files, the helper displays: All three comparison products are ready. An HTTP 401 response means the token needs to be refreshed using the authentication procedure in the prerequisites. The same download command can then be repeated; complete extracted products are reused.

The download is complete when this input directory contains three extracted .SAFE folders, one for each acquisition date, and comparison-products.json. The three ZIP archives are also kept there:

Downloads/aral-data/
├── comparison-products.json
├── S2A_MSIL2A_20230526T065621_N0510_R063_T40TGS_20240907T073145.SAFE/
├── S2A_MSIL2A_20230605T065621_N0510_R063_T40TGS_20240911T135334.SAFE/
└── S2A_MSIL2A_20230705T065631_N0510_R063_T40TGS_20240912T130104.SAFE/

The folder name is part of the input configuration: the notebook expects aral-data directly inside the current user’s Downloads folder. Each full .SAFE directory name must match the path in comparison-products.json, including its spelling, case and suffix. A ZIP archive cannot be used as a processing directory. The files MTD_MSIL2A.xml and GRANULE must be directly inside each extracted .SAFE folder.

For a different location, the helper’s –data-dir argument and the notebook’s manifest path must name that same directory. Rerunning the helper after moving the products updates the absolute paths recorded in the manifest. The input folder aral-data and output folder aral_ndwi have distinct purposes and names.

Obtaining NDWI

Convert B03 and B08 digital numbers to surface reflectance before calculating NDWI. The per-band offsets introduced with Processing Baseline 04.00 do not cancel in the denominator of the index. Floating-point arithmetic also avoids integer overflow when adding or subtracting band values.

from pathlib import Path
import xml.etree.ElementTree as ET

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.mask import mask
from rasterio.vrt import WarpedVRT
from rasterio.warp import transform_geom

Functions for reading, calculating and saving raster data

The functions below read the quantification value and B03/B08 offsets from MTD_MSIL2A.xml. They crop each granule to the bounding box and align the 20 m SCL layer with the 10 m bands using nearest-neighbour resampling. Only SCL classes 4, 5 and 6 (vegetation, bare soil and water) are retained. No-data pixels, negative reflectance and nonpositive denominators are excluded. Missing calibration metadata stops processing rather than silently using incorrect values.

def read_scaling(product_dir):
    metadata = product_dir / "MTD_MSIL2A.xml"
    root = ET.parse(metadata).getroot()
    elements = list(root.iter())

    def value(tag):
        return next(e.text for e in elements if e.tag.split("}")[-1] == tag)

    quantification = float(value("BOA_QUANTIFICATION_VALUE"))
    baseline = float(value("PROCESSING_BASELINE"))
    offsets = {
        e.attrib["band_id"]: float(e.text)
        for e in elements if e.tag.split("}")[-1] == "BOA_ADD_OFFSET"
    }
    # Sentinel-2 metadata numbers B03 and B08 as band IDs 2 and 7.
    if baseline < 4.0 and not offsets:
        offsets = {"2": 0.0, "7": 0.0}
    if not np.isfinite(quantification) or quantification <= 0:
        raise ValueError(f"Invalid reflectance quantification in {metadata}")
    if any(key not in offsets or not np.isfinite(offsets[key]) for key in ("2", "7")):
        raise ValueError(f"Missing or invalid B03/B08 offsets in {metadata}")
    return quantification, offsets


def find_band(granule, band, resolution):
    matches = sorted(granule.glob(f"IMG_DATA/R{resolution}m/*_{band}_{resolution}m.jp2"))
    if len(matches) != 1:
        raise ValueError(f"Expected one {band} file at {resolution} m in {granule}")
    return matches[0]


def process_granule(product_dir, granule, bbox, output_dir):
    quantification, offsets = read_scaling(product_dir)
    minx, miny, maxx, maxy = bbox
    geometry = {
        "type": "Polygon",
        "coordinates": [[[minx, miny], [maxx, miny], [maxx, maxy],
                         [minx, maxy], [minx, miny]]],
    }
    with (
        rasterio.open(find_band(granule, "B03", 10)) as green_src,
        rasterio.open(find_band(granule, "B08", 10)) as nir_src,
        rasterio.open(find_band(granule, "SCL", 20)) as scl_src,
    ):
        if green_src.crs is None or scl_src.crs is None:
            raise ValueError(f"Missing raster CRS in {granule}")
        if (green_src.crs != nir_src.crs or green_src.transform != nir_src.transform
                or green_src.shape != nir_src.shape):
            raise ValueError(f"B03 and B08 grids differ in {granule}")
        shape = transform_geom("EPSG:4326", green_src.crs, geometry)
        try:
            green_dn, transform = mask(green_src, [shape], crop=True, filled=False, indexes=1)
        except ValueError as exc:
            if "do not overlap" in str(exc):
                return None
            raise
        nir_dn, nir_transform = mask(nir_src, [shape], crop=True, filled=False, indexes=1)
        if transform != nir_transform or green_dn.shape != nir_dn.shape:
            raise ValueError(f"Cropped B03 and B08 grids differ in {granule}")
        # Resample the categorical 20 m SCL grid onto the 10 m band grid.
        with WarpedVRT(scl_src, crs=green_src.crs, transform=green_src.transform,
                       width=green_src.width, height=green_src.height,
                       resampling=Resampling.nearest) as scl_vrt:
            scl, _ = mask(scl_vrt, [shape], crop=True, filled=False, indexes=1)

        green = (green_dn.data.astype("float32") + offsets["2"]) / quantification
        nir = (nir_dn.data.astype("float32") + offsets["7"]) / quantification
        denominator = green + nir
        valid = (
            ~np.ma.getmaskarray(green_dn) & ~np.ma.getmaskarray(nir_dn)
            & ~np.ma.getmaskarray(scl)
            & (green_dn.data != 0) & (nir_dn.data != 0)
            & np.isin(scl.data, [4, 5, 6])  # vegetation, bare soil, water
            & np.isfinite(green) & np.isfinite(nir)
            & (green >= 0) & (nir >= 0) & (denominator > 0)
        )
        if not valid.any():
            return None
        ndwi = np.full(green.shape, -9999.0, dtype="float32")
        np.divide(green - nir, denominator, out=ndwi, where=valid)
        output = output_dir / f"{product_dir.stem}_{granule.name}_NDWI.tif"
        with rasterio.open(output, "w", driver="GTiff", height=ndwi.shape[0],
                           width=ndwi.shape[1], count=1, dtype="float32",
                           crs=green_src.crs, transform=transform,
                           nodata=-9999.0, compress="deflate") as dst:
            dst.write(ndwi, 1)
    return output

Select local products and compute NDWI rasters

The manifest contains the three dates and the absolute paths of the extracted local products. This cell sets dataset to those paths in date order. Its output should show 2023-05-26, 2023-06-05 and 2023-07-05.

import json
from pathlib import Path

manifest = Path.home() / "Downloads/aral-data/comparison-products.json"
comparison_products = json.loads(manifest.read_text())

dataset = [item["path"] for item in comparison_products]

for item in comparison_products:
    print(item["date"], item["name"])

The selection cell runs immediately before the processing loop. A dataset containing only one product produces only one raster, which is insufficient for the three-panel comparison. The parameter, import and function-definition cells must have run in the current kernel session.

The loop creates aral_ndwi in the notebook’s current working directory and prints three GeoTIFF paths when processing finishes. Each output retains its granule’s native projected grid and uses -9999 as its no-data value. Granules outside the bounding box or without valid pixels are skipped; overlapping tiles are not mosaicked.

output_dir = Path.cwd() / "aral_ndwi"
output_dir.mkdir(parents=True, exist_ok=True)
outputs = []

for item in dict.fromkeys(dataset):
    product_dir = Path(item)
    if not product_dir.is_dir():
        raise FileNotFoundError(
            f"Product is not accessible: {product_dir}. "
            "Check the extracted .SAFE folder and its path in comparison-products.json."
        )
    granules = sorted(p for p in (product_dir / "GRANULE").iterdir() if p.is_dir())
    if not granules:
        raise ValueError(f"No granules found in {product_dir}")
    for granule in granules:
        output = process_granule(product_dir, granule, bbox, output_dir)
        if output is not None:
            outputs.append(output)
            print(output)

if not outputs:
    raise RuntimeError("No valid pixels remain inside the bounding box after masking.")

Display the first generated raster with a fixed NDWI colour scale. Masked pixels are transparent. Select another entry in outputs to inspect a different raster.

import matplotlib.pyplot as plt
from rasterio.plot import plotting_extent

with rasterio.open(outputs[0]) as src:
    ndwi = src.read(1, masked=True)
    fig, ax = plt.subplots(figsize=(10, 6))
    image = ax.imshow(ndwi, extent=plotting_extent(src), cmap="RdYlBu",
                      vmin=-1, vmax=1)
    ax.set_title(outputs[0].name)
    ax.set_xlabel(f"Easting ({src.crs})")
    ax.set_ylabel(f"Northing ({src.crs})")
    fig.colorbar(image, ax=ax, label="NDWI")
    plt.show()

Compare and interpret the results

The comparison uses the three outputs for 26 May, 5 June and 5 July 2023, in the same order as the manifest. All three selected products cover tile T40TGS from relative orbit R063.

The plotting code checks that the rasters have matching coordinate systems, pixel grids and dimensions. Reduced arrays are read for display; the GeoTIFFs remain at full resolution. The panels share the same map extent and NDWI scale. A common mask leaves a displayed location white if any of the three dates lacks a valid observation there, so white areas do not indicate dry land.

import numpy as np
import rasterio
import matplotlib.pyplot as plt
from rasterio.enums import Resampling

if len(outputs) != 3:
    raise ValueError("Run section 3.3 with the three comparison products first.")

images = []
reference_grid = None

for path in outputs:
    with rasterio.open(path) as src:
        grid = (src.crs, src.transform, src.shape)
        if reference_grid is None:
            reference_grid = grid
            extent = (
                src.bounds.left, src.bounds.right,
                src.bounds.bottom, src.bounds.top,
            )
            # Read a smaller version for display; the GeoTIFF stays unchanged.
            height = min(1200, src.height)
            width = round(src.width * height / src.height)
        elif grid != reference_grid:
            raise ValueError("The rasters must use the same grid.")

        images.append(src.read(
            1, masked=True, out_shape=(height, width),
            resampling=Resampling.nearest,
        ))

# Display only locations with valid observations on all three dates.
common_mask = np.logical_or.reduce([
    np.ma.getmaskarray(image) for image in images
])

fig, axes = plt.subplots(
    1, 3, figsize=(16, 6),
    sharex=True, sharey=True, layout="constrained",
)

for ax, image, item in zip(axes, images, comparison_products):
    displayed = np.ma.array(image.data, mask=common_mask)
    plot = ax.imshow(
        displayed, extent=extent, cmap="RdYlBu",
        vmin=-1, vmax=1, interpolation="nearest",
    )
    ax.set_title(item["date"])
    ax.set_xlabel("Easting (m)")

axes[0].set_ylabel("Northing (m)")
fig.colorbar(plot, ax=axes, label="NDWI", shrink=0.8)
plt.show()

Comparison generated from the downloaded products using the reflectance conversion and quality mask above.

The three-raster check stops plotting if outputs has a different number of entries. The manifest-selection cell and processing loop must then run again before plotting. A kernel restart also requires the parameters, imports and function definitions to be rerun.

Blue shades represent higher NDWI values, while yellow and orange shades represent lower values. Differences between dates describe changes in green and near-infrared reflectance contrast at locations with valid observations. The colour scale supports a visual comparison; it does not establish water depth, volume or the cause of a change.

These are three individual acquisitions, not daily records or monthly composites. The common mask applies to the reduced display arrays. Quantitative estimates of water extent would require a justified classification method and a common valid-data mask on the full-resolution grids.

The result of this workflow is three georeferenced NDWI rasters and a comparison of the same area on a common colour scale, with locations lacking valid observations excluded from all three panels. You now have a repeatable way to inspect apparent changes in water coverage near the river mouth and flag locations for closer investigation. The saved GeoTIFFs are also available for further mapping in a GIS. Together, these outputs provide a starting point for planning field checks or comparing the satellite observations with river-flow and water-level records.

What to do next

In this article, we downloaded three Sentinel-2 products as ZIP archives and processed the extracted files locally to compare acquisitions from different dates. Other ways to access satellite data include an EODATA filesystem mount or an S3 client. The guides below explain how to set up these alternatives.

s3fs provides filesystem access to the EODATA repository, usually through /eodata. Applications read satellite files through that mount without first downloading complete ZIP archives. On a Creodias cloud VM, How to mount eodata using s3fs in Linux on Creodias describes the setup, including the required EODATA network access and S3 credentials.

Other clients, including boto3 and s3cmd, can read or copy satellite data through the S3 API. The available methods and credential setup are described in How to access EODATA using S3 protocol on Creodias.

S3 access from a local workstation requires the external EODATA service and its endpoint and credentials, described in How to access EODATA from your own infrastructure on Creodias. The access token used for ZIP downloads does not replace S3 access keys.

Interested to see additional examples of studying nature via satellite images?

Urban sprawl:

Monitoring Urban Sprawl with Sentinel-1 SAR Data Using Pixel Value and Polarization Thresholding on Creodias

Floods:

Analyzing and monitoring floods using Python and Sentinel-2 satellite imagery on Creodias