Cloud-Native Intertidal Vegetation Monitoring with Sentinel-2 and EOPF Zarr

By: @fer-at-sea and @pglez82

Introduction

Intertidal zones β€” the strip of coastline exposed at low tide and submerged at high tide β€” are among the most productive and ecologically significant habitats on Earth. Despite covering less than 0.5 % of the global ocean surface, they provide key ecosystem services:

  • Carbon sequestration: Salt marshes and seagrass meadows are β€œblue carbon” ecosystems that can sequester carbon at rates several times higher than terrestrial forests and store much of it in long-term sediments.
  • Coastal protection: Intertidal vegetation dissipates wave energy, reducing flood risk for coastal communities in a warming climate.
  • Biodiversity hotspots: These zones support invertebrates, migratory shorebirds, juvenile fish, and endemic plant communities that exist nowhere else.
  • Climate sensitivity: Because intertidal vegetation is simultaneously exposed to the atmosphere (temperature, drought) and the sea (salinity, sea-level change), it acts as an early-warning indicator of climate-driven coastal stress.

Study Area

The Villaviciosa estuary (Asturias, northern Spain) is a Ramsar-designated wetland of international importance hosting extensive intertidal mudflats and saltmarshes. Regular and systematic monitoring at fine spatial scales can provide valuable information on habitat dynamics but is often difficult to sustain using field surveys alone. Satellite remote sensing, particularly within cloud-native processing frameworks, offers a complementary approach that can facilitate basin-scale and time-series monitoring.

Key characteristics relevant to our remote-sensing approach:

Feature Detail
Tidal regime Semi-diurnal; tidal range ~4 m at spring tides
Exposures Extensive low-gradient mudflats exposed for 6–8 h per tidal cycle
Designations Ramsar site, Natura 2000 SAC/SPA, Asturias Nature Reserve
Sentinel-2 tiles Primarily 29TQJ / 30TUP (border zone between UTM zones 29N and 30N)

A semi-diurnal tidal regime with ~4 m range means that at any given Sentinel-2 overpass (roughly every 5 days), the estuary may be near fully exposed, near fully submerged, or anywhere in between. This tidal variability is the key driver of the water-frequency approach: rather than attempting to record a single synoptic state, we accumulate many snapshots across the full tidal cycle and use the fraction of time a pixel is wet as a proxy for intertidal zone extent.

What we will learn

  • πŸ›°οΈ How to discover Sentinel-2 products via EOPF STAC and open them as Zarr with xarray.
  • 🌊 How to build a time series over a coastal area using lazy loading and partial reads (cloud-native workflow).
  • 🌿 How to derive simple spectral indicators to support intertidal vegetation monitoring and visualise the results.

Prerequisites

Workflow Overview

We build a complete cloud-native workflow in six stages:

  1. Data discovery β€” query the EOPF STAC catalogue for Sentinel-2 L2A tiles covering the study area using pystac_client, with no data download required.
  2. Lazy cloud access β€” open each tile as a Dask-backed xarray dataset directly from object storage via the eopf-zarr engine and clip it to the area of interest.
  3. Water frequency mapping β€” stack the Scene Classification Layer (SCL) across all acquisitions to compute a per-pixel water frequency, exploiting the tidal cycle as a natural inundation signal.
  4. Intertidal zone delineation β€” threshold the water frequency map to isolate intermittently flooded pixels and export the result as a GeoTIFF and GeoPackage vector file.
  5. Vegetation monitoring β€” for each cloud-free acquisition, compute the Normalised Difference Vegetation Index (NDVI) from the red (B04) and near-infrared (B08) bands, restrict the analysis to exposed intertidal pixels (SCL 4/5), and track the fraction with NDVI > 0.2 over time. Normalising by the number of exposed pixels rather than total intertidal area removes the tidal confounding effect, giving a tide-independent vegetation cover metric.
  6. Temporal aggregation β€” summarise the NDVI time series by calendar month and by season (Winter / Spring / Summer / Autumn) to reveal phenological patterns in intertidal vegetation cover.

A note on data availability The EOPF Zarr Sample Service STAC catalogue maintains a rolling archive of ~1 year for Europe. This has two important implications for this workflow:

  • Intertidal mask quality improves with more images. A longer time series captures a wider range of tidal states β€” including extreme high and low tides β€” which helps the water-frequency approach converge on the true intertidal boundary. With only a year of data the mask may be incomplete or biased towards the tidal conditions that happened to be observed.
  • NDVI temporal trend becomes more meaningful as the archive grows. Seasonal vegetation cycles (e.g. growth of seagrass or macroalgal blooms) can only be detected once multi-month records are available. As the rolling window advances, re-running this notebook will automatically incorporate new acquisitions.

Import libraries

# Standard library
from concurrent.futures import ThreadPoolExecutor, as_completed

# Scientific computing
import numpy as np
import pandas as pd

# Geospatial 
import rasterio
import rasterio.enums
import rasterio.transform
from rasterio.features import shapes
from rasterio.warp import reproject as rio_reproject, Resampling, calculate_default_transform
import rioxarray  # noqa: F401  – required to activate the .rio accessor on xarray objects
import xarray as xr
import geopandas as gpd
from shapely.geometry import box, shape
from pyproj import Transformer
from pyproj import CRS

# Visualisation
import matplotlib.pyplot as plt
import contextily as ctx
import warnings

# Target the ZarrUserWarning explicitly by its message pattern
warnings.filterwarnings("ignore", message=".*Both zarr.json.*and .zgroup.*")
warnings.filterwarnings("ignore", message=".*Analysis mode not implemented.*")

# EO data access
import pystac_client

Area of interest

site = "Villaviciosa"
bbox = [-5.44, 43.48, -5.36, 43.54]
crs_wgs84 = "EPSG:4326"
crs_web_mercator = "EPSG:3857"


# -----------------------------
# STUDY AREA MAP
# -----------------------------
aoi_wm = gpd.GeoDataFrame(geometry=[box(*bbox)], crs=crs_wgs84).to_crs(crs_web_mercator)

fig, (ax_ov, ax_det) = plt.subplots(1, 2, figsize=(14, 6))

# ---- Left ----
spain_wm = gpd.GeoDataFrame(
    geometry=[box(-9.5, 36.0, 4.5, 44.2)], crs=crs_wgs84
).to_crs(crs_web_mercator)
b = spain_wm.total_bounds
ax_ov.set_xlim(b[0], b[2])
ax_ov.set_ylim(b[1], b[3])
ctx.add_basemap(ax_ov, crs=crs_web_mercator, source=ctx.providers.CartoDB.Positron, zoom=6)
aoi_wm.boundary.plot(ax=ax_ov, color="red", linewidth=1.5)
centroid = aoi_wm.geometry.iloc[0].centroid
ax_ov.plot(centroid.x, centroid.y, marker="*", color="red", markersize=14, zorder=5)
ax_ov.set_axis_off()

# ---- Right ----
BUFFER = 2500   
b = aoi_wm.total_bounds
ax_det.set_xlim(b[0] - BUFFER, b[2] + BUFFER)
ax_det.set_ylim(b[1] - BUFFER, b[3] + BUFFER)
ctx.add_basemap(ax_det, crs=crs_web_mercator, source=ctx.providers.CartoDB.Positron, zoom=13)
aoi_wm.boundary.plot(ax=ax_det, color="red", linewidth=2, label="AOI")
ax_det.legend(loc="lower right", fontsize=9)
ax_det.set_axis_off()

fig.suptitle(f"Study Area: {site}", fontsize=13, fontweight="bold")
plt.tight_layout()
plt.show()

Fetching S2 L2A zarr data

Load into xarray

Why eopf-zarr + Dask?

Traditional EO workflows require downloading full granules (often hundreds of MB per scene) before any analysis. The eopf-zarr xarray engine replaces that pattern with streaming access: data is read in small HTTP range-requests from the object store only when a computation actually requires it. Wrapping the dataset with Dask (chunks={}) defers even those reads until explicitly triggered β€” the cell below executes relatively fast, because it only fetches metadata and coordinates.

Key concepts: - xr.open_dataset(..., engine="eopf-zarr", chunks={}) β€” opens the store lazily; all variable arrays are Dask graphs, not in-memory arrays. - .sel(x=..., y=...) β€” label-based slicing restricts subsequent reads to the AOI bounding box; pixels outside the window are never fetched. - CRS alignment β€” the study area straddles the UTM 29N / 30N tile boundary, so we reproject the WGS-84 bounding box into each tile’s native CRS before slicing.

For each STAC item we open the corresponding Zarr store with the eopf-zarr engine β€” no file download required. We then clip each dataset to the AOI using xarray label-based indexing and skip any item that has no overlap or is missing a CRS. Dask chunks (chunks={}) ensure that data is read lazily and only the pixels within the AOI are actually loaded into memory.

Sequential vs. parallel loading: The notebook provides two alternative cells β€” a simple sequential loop (easy to follow) and a concurrent version using ThreadPoolExecutor. The parallel version issues all Zarr metadata requests simultaneously, dramatically reducing wall-clock time when the archive contains many items. The outputs of both cells are identical; only one needs to be run.

# datasets, included_items, skipped, crs_list = [], [], [], []

# for item in items:
#     url = item.assets["product"].href

#     # --- Try to open dataset ---
#     try:
#         ds = xr.open_dataset(url, engine="eopf-zarr", chunks={})
#     except Exception as e:
#         print(f"\n-{e}")
#         skipped.append(url)
#         continue

#     # --- Get CRS ---
#     crs = ds.rio.crs
#     print(f"\n-{url.split('/')[-1]}")
#     print(f"   CRS: {crs}")

#     if crs is None:
#         print("   No CRS found, skipping")
#         skipped.append(url)
#         continue

#     # --- Reproject WGS84 bbox to dataset CRS ---
#     transformer = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
#     minx, miny = transformer.transform(bbox[0], bbox[1])
#     maxx, maxy = transformer.transform(bbox[2], bbox[3])

#     # --- Handle Y axis direction ---
#     y_slice = slice(maxy, miny) if ds.y[0] > ds.y[-1] else slice(miny, maxy)

#     # --- Clip ---
#     ds_clipped = ds.sel(x=slice(minx, maxx), y=y_slice)

#     if ds_clipped.sizes["x"] == 0 or ds_clipped.sizes["y"] == 0:
#         print("   No overlap with bbox after reprojection")
#         skipped.append(url)
#         continue

#     print(f"   Clipped size: {dict(ds_clipped.sizes)}")

#     # --- Store dataset and its STAC item together ---
#     datasets.append(ds_clipped)
#     included_items.append(item)
#     crs_list.append(crs)

# print(f"\nβœ“ Included: {len(included_items)} | βœ— Skipped: {len(skipped)}")

The following cell is a Parallel version of the cell above, aiming to run all Zarr items concurrently. The results are sorted by acquisition date so downstream cells are unaffected.

MAX_WORKERS = 8   # concurrent metadata fetches; tune to network conditions

# Transformer cache: avoids rebuilding the same CRS projection per thread.
# NW Spain tiles appear in both EPSG:32629 (Zone 29N) and EPSG:32630 (Zone 30N).
_transformer_cache: dict = {}

def _open_and_clip(item):
    """Open the Zarr store, isolate ONLY required variables, and clip immediately."""
    url = item.assets["product"].href
    tile_name = url.split('/')[-1]
    msgs = []

    # ---- FIRST FILTER: OPENING THE DATA ----
    try:
        dt = xr.open_datatree(url, engine="eopf-zarr", chunks={})
        
        # Isolate variables lazily
        ds_ref = dt["measurements/reflectance/r20m"]['b04'].to_dataset()
        ds_mask = dt["conditions/mask/l2a_classification/r20m"]['scl'].to_dataset()
        ds_10 = dt["measurements/reflectance/r10m"]['b08'].to_dataset()
        
    except Exception as e:
        # EXPLICIT LOGGING: Show it failed the first filter immediately
        print(f"❌ SKIP {tile_name} β€” Failed First Filter (Could not open or locate bands). Error: {e}")
        msgs.append(f"   ERROR opening DataTree: Skip")
        return None, item, None, msgs  # Force an exit right here!
    
    # ---- SECOND FILTER: CRS CHECK ----
    raw_crs = dt.attrs.get('stac_discovery', {}).get('properties', {}).get('proj:code')
    if not raw_crs:
        print(f"❌ SKIP {tile_name} β€” Failed Second Filter (Missing 'proj:code' in metadata)")
        msgs.append("   No CRS string found in metadata, skipping")
        return None, item, None, msgs

    try:
        crs = CRS.from_user_input(raw_crs)
        
        # Inject CRS configuration
        for d in [ds_ref, ds_mask, ds_10]:
            d.rio.set_spatial_dims(x_dim="x", y_dim="y", inplace=True)
            d.rio.write_crs(crs, inplace=True)
        
        # Coordinate math
        crs_key = str(crs)
        if crs_key not in _transformer_cache:
            _transformer_cache[crs_key] = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
        transformer = _transformer_cache[crs_key]

        minx, miny = transformer.transform(bbox[0], bbox[1])
        maxx, maxy = transformer.transform(bbox[2], bbox[3])
        
        y_slice = slice(maxy, miny) if ds_ref.y[0] > ds_ref.y[-1] else slice(miny, maxy)
        
        # ---- THIRD FILTER: BOUNDING BOX MATCH ----
        ds_ref_clip = ds_ref.sel(x=slice(minx, maxx), y=y_slice)
        ds_mask_clip = ds_mask.sel(x=slice(minx, maxx), y=y_slice)
        ds_10_clip = ds_10.sel(x=slice(minx, maxx), y=y_slice)

        if ds_ref_clip.sizes["x"] == 0 or ds_ref_clip.sizes["y"] == 0:
            print(f"❌ SKIP {tile_name} β€” Failed Third Filter (No overlap with bbox)")
            return None, item, None, msgs

        # Resample and merge tiny clips
        mask_10m_res = ds_10_clip.rio.reproject_match(ds_ref_clip, resampling=rasterio.enums.Resampling.nearest)
        mask_10m_res = mask_10m_res.rename({v: f"{v}_10m" for v in mask_10m_res.data_vars})
                
        ds_final = xr.merge([ds_ref_clip, ds_mask_clip, mask_10m_res], compat="override")
        
        # SUCCESS LOG
        print(f"βœ… INCLUDE {tile_name} β€” Successfully parsed and clipped")
        
    except Exception as err:
        print(f"❌ SKIP {tile_name} β€” Failed Processing Filter. Error: {err}")
        msgs.append(f"   Processing Error")
        return None, item, None, msgs  # Ensure we exit on error!

    return ds_final, item, crs, msgs

# Set variables to store results 

datasets, included_items, skipped, crs_list = [], [], [], []

with ThreadPoolExecutor(max_workers=MAX_WORKERS) as pool:
    future_to_item = {pool.submit(_open_and_clip, item): item for item in items}
    raw_results = []
    for future in as_completed(future_to_item):
        ds_clipped, item, crs, msgs = future.result()
        if ds_clipped is None and any("ERROR" in m for m in msgs):
                    for m in msgs:
                        print(m)
        raw_results.append((item.datetime, ds_clipped, item, crs))

# Sort chronologically β€” as_completed returns in arrival order, not submission order
raw_results.sort(key=lambda r: r[0])

for _, ds_clipped, item, crs in raw_results:
    if ds_clipped is not None:
        datasets.append(ds_clipped)
        included_items.append(item)
        crs_list.append(crs)
    else:
        skipped.append(item.assets["product"].href)

print(f"\nβœ“ Included: {len(included_items)} | βœ— Skipped: {len(skipped)}")
βœ… INCLUDE S2A_MSIL2A_20260321T112131_N0512_R037_T29TQJ_20260321T200016.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260321T112131_N0512_R037_T30TUP_20260321T200016.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260331T112131_N0512_R037_T30TUP_20260331T200509.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260321T110649_N0512_R137_T30TUP_20260321T134453.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260324T112109_N0512_R037_T30TUP_20260324T150715.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260329T112111_N0512_R037_T29TQJ_20260329T144222.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260329T112111_N0512_R037_T30TUP_20260329T144222.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260324T112109_N0512_R037_T29TQJ_20260324T150715.zarr β€” Successfully parsed and clipped
❌ SKIP S2A_MSIL2A_20260311T112131_N0512_R037_T29TQJ_20260311T201411.zarr β€” Failed Second Filter (Missing 'proj:code' in metadata)
βœ… INCLUDE S2B_MSIL2A_20260321T110649_N0512_R137_T29TQJ_20260321T134453.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260319T112121_N0512_R037_T30TUP_20260319T164912.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260319T112121_N0512_R037_T29TQJ_20260319T164912.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260318T110731_N0512_R137_T30TUP_20260318T180922.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260316T110821_N0512_R137_T29TQJ_20260316T152432.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260318T110731_N0512_R137_T29TQJ_20260318T180922.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260316T110821_N0512_R137_T30TUP_20260316T152432.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260308T110841_N0512_R137_T30TUP_20260308T175910.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260304T112109_N0512_R037_T29TQJ_20260304T165437.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260308T110841_N0512_R137_T29TQJ_20260308T175910.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260301T110909_N0512_R137_T30TUP_20260301T150305.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260304T112109_N0512_R037_T30TUP_20260304T165437.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260301T110909_N0512_R137_T29TQJ_20260301T150305.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260224T111051_N0512_R137_T30TUP_20260224T145311.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260224T111051_N0512_R137_T29TQJ_20260224T145311.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260222T112109_N0512_R037_T29TQJ_20260222T165721.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260222T112109_N0512_R037_T30TUP_20260222T165721.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260219T111621_N0512_R137_T30TUP_20260219T163114.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260219T111009_N0512_R137_T30TUP_20260219T134643.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20260219T111621_N0512_R137_T29TQJ_20260219T163114.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260219T111009_N0512_R137_T29TQJ_20260219T134643.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260217T112131_N0512_R037_T30TUP_20260217T150410.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260214T111151_N0512_R137_T30TUP_20260214T145911.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260214T111151_N0512_R137_T29TQJ_20260214T145911.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260212T112109_N0512_R037_T30TUP_20260212T151217.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260212T112109_N0512_R037_T29TQJ_20260212T151217.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260123T112249_N0511_R037_T30TUP_20260123T133836.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20260118T112411_N0511_R037_T29TQJ_20260118T130309.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20260110T111339_N0511_R137_T30TUP_20260110T131623.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251231T111359_N0511_R137_T30TUP_20251231T134320.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251231T111359_N0511_R137_T29TQJ_20251231T134320.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251229T112511_N0511_R037_T30TUP_20251229T125809.zarr β€” Successfully parsed and clipped
❌ SKIP S2B_MSIL2A_20251211T111349_N0511_R137_T30TUP_20251212T112638.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202512-s02msil2a-eu/11/products/cpm_v262/S2B_MSIL2A_20251211T111349_N0511_R137_T30TUP_20251212T112638.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2B_MSIL2A_20251211T111349_N0511_R137_T30TUP_20251211T120627.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202512-s02msil2a-eu/11/products/cpm_v262/S2B_MSIL2A_20251211T111349_N0511_R137_T30TUP_20251211T120627.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2B_MSIL2A_20251211T111349_N0511_R137_T29TQJ_20251212T112638.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202512-s02msil2a-eu/11/products/cpm_v262/S2B_MSIL2A_20251211T111349_N0511_R137_T29TQJ_20251212T112638.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2B_MSIL2A_20251211T111349_N0511_R137_T29TQJ_20251211T120627.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202512-s02msil2a-eu/11/products/cpm_v262/S2B_MSIL2A_20251211T111349_N0511_R137_T29TQJ_20251211T120627.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2A_MSIL2A_20251228T111521_N0511_R137_T29TQJ_20251228T131811.zarr β€” Successfully parsed and clipped
❌ SKIP S2A_MSIL2A_20251208T111501_N0511_R137_T29TQJ_20251208T131613.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202512-s02msil2a-eu/08/products/cpm_v262/S2A_MSIL2A_20251208T111501_N0511_R137_T29TQJ_20251208T131613.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2B_MSIL2A_20251214T112359_N0511_R037_T30TUP_20251214T121801.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251226T111501_N0511_R137_T30TUP_20251226T124609.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251221T112521_N0511_R037_T30TUP_20251221T151311.zarr β€” Successfully parsed and clipped
❌ SKIP S2C_MSIL2A_20251126T111421_N0511_R137_T30TUP_20251126T151514.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/26/products/cpm_v262/S2C_MSIL2A_20251126T111421_N0511_R137_T30TUP_20251126T151514.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2B_MSIL2A_20251221T111359_N0511_R137_T29TQJ_20251221T134743.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251214T112359_N0511_R037_T29TQJ_20251214T121801.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251211T112511_N0511_R037_T30TUP_20251211T151014.zarr β€” Successfully parsed and clipped
❌ SKIP S2A_MSIL2A_20251121T112421_N0511_R037_T30TUP_20251121T151713.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/21/products/cpm_v262/S2A_MSIL2A_20251121T112421_N0511_R037_T30TUP_20251121T151713.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2A_MSIL2A_20251121T112421_N0511_R037_T29TQJ_20251121T151713.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/21/products/cpm_v262/S2A_MSIL2A_20251121T112421_N0511_R037_T29TQJ_20251121T151713.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2B_MSIL2A_20251121T111259_N0511_R137_T30TUP_20251121T120801.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/21/products/cpm_v262/S2B_MSIL2A_20251121T111259_N0511_R137_T30TUP_20251121T120801.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2B_MSIL2A_20251121T111259_N0511_R137_T29TQJ_20251121T120801.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/21/products/cpm_v262/S2B_MSIL2A_20251121T111259_N0511_R137_T29TQJ_20251121T120801.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2C_MSIL2A_20251206T111451_N0511_R137_T29TQJ_20251206T124813.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251206T111451_N0511_R137_T30TUP_20251206T124813.zarr β€” Successfully parsed and clipped
❌ SKIP S2C_MSIL2A_20251119T112401_N0511_R037_T29TQJ_20251119T130615.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/19/products/cpm_v262/S2C_MSIL2A_20251119T112401_N0511_R037_T29TQJ_20251119T130615.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2A_MSIL2A_20251128T111431_N0511_R137_T29TQJ_20251128T132409.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251128T111431_N0511_R137_T30TUP_20251128T132409.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251126T111421_N0511_R137_T29TQJ_20251126T151514.zarr β€” Successfully parsed and clipped
❌ SKIP S2A_MSIL2A_20251108T111311_N0511_R137_T29TQJ_20251108T150910.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/08/products/cpm_v262/S2A_MSIL2A_20251108T111311_N0511_R137_T29TQJ_20251108T150910.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2B_MSIL2A_20251101T111119_N0511_R137_T29TQJ_20251101T135730.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202511-s02msil2a-eu/01/products/cpm_v262/S2B_MSIL2A_20251101T111119_N0511_R137_T29TQJ_20251101T135730.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2B_MSIL2A_20251124T112309_N0511_R037_T30TUP_20251124T133211.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251124T112309_N0511_R037_T29TQJ_20251124T133211.zarr β€” Successfully parsed and clipped
❌ SKIP S2A_MSIL2A_20251019T111111_N0511_R137_T30TUP_20251019T151521.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202510-s02msil2a-eu/19/products/cpm_v256/S2A_MSIL2A_20251019T111111_N0511_R137_T30TUP_20251019T151521.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2C_MSIL2A_20251119T112401_N0511_R037_T30TUP_20251119T130615.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251116T111341_N0511_R137_T30TUP_20251116T124312.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251108T111311_N0511_R137_T30TUP_20251108T150910.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251109T112321_N0511_R037_T30TUP_20251109T130709.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251109T112321_N0511_R037_T29TQJ_20251109T130709.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251027T111201_N0511_R137_T30TUP_20251027T145211.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251027T111201_N0511_R137_T29TQJ_20251027T145211.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251019T111111_N0511_R137_T29TQJ_20251019T151521.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251017T111101_N0511_R137_T30TUP_20251017T151300.zarr β€” Successfully parsed and clipped
❌ SKIP S2A_MSIL2A_20251009T111001_N0511_R137_T30TUP_20251009T145813.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202510-s02msil2a-eu/09/products/cpm_v256/S2A_MSIL2A_20251009T111001_N0511_R137_T30TUP_20251009T145813.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2C_MSIL2A_20251017T111101_N0511_R137_T29TQJ_20251017T151300.zarr β€” Successfully parsed and clipped
❌ SKIP S2C_MSIL2A_20251007T110951_N0511_R137_T30TUP_20251007T145121.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202510-s02msil2a-eu/07/products/cpm_v256/S2C_MSIL2A_20251007T110951_N0511_R137_T30TUP_20251007T145121.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2B_MSIL2A_20251015T112109_N0511_R037_T29TQJ_20251015T133117.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251012T112131_N0511_R037_T30TUP_20251012T152415.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251012T112131_N0511_R037_T29TQJ_20251012T152415.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251012T110919_N0511_R137_T30TUP_20251012T133527.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251012T110919_N0511_R137_T29TQJ_20251012T133527.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20251010T112131_N0511_R037_T29TQJ_20251010T151313.zarr β€” Successfully parsed and clipped
❌ SKIP S2B_MSIL2A_20250925T112109_N0511_R037_T30TUP_20250925T151904.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202509-s02msil2a-eu/25/products/cpm_v256/S2B_MSIL2A_20250925T112109_N0511_R037_T30TUP_20250925T151904.zarr' at path ''
   ERROR opening DataTree: Skip
❌ SKIP S2B_MSIL2A_20250925T112109_N0511_R037_T29TQJ_20250925T151904.zarr β€” Failed First Filter (Could not open or locate bands). Error: No group found in store 'https://objects.eodc.eu:443/e05ab01a9d56408d82ac32d69a5aae2a:202509-s02msil2a-eu/25/products/cpm_v256/S2B_MSIL2A_20250925T112109_N0511_R037_T29TQJ_20250925T151904.zarr' at path ''
   ERROR opening DataTree: Skip
βœ… INCLUDE S2A_MSIL2A_20251009T111001_N0511_R137_T29TQJ_20251009T145813.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251002T112131_N0511_R037_T30TUP_20251002T152717.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20251002T112131_N0511_R037_T29TQJ_20251002T152717.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251002T110809_N0511_R137_T29TQJ_20251002T135509.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20251002T110809_N0511_R137_T30TUP_20251002T135509.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250929T110851_N0511_R137_T30TUP_20250929T151314.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250929T110851_N0511_R137_T29TQJ_20250929T151314.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250919T110741_N0511_R137_T30TUP_20250919T152715.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250922T112131_N0511_R037_T29TQJ_20250922T160420.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250919T110741_N0511_R137_T29TQJ_20250919T152715.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250917T110731_N0511_R137_T29TQJ_20250917T163216.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250917T110731_N0511_R137_T30TUP_20250917T163216.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250907T110631_N0511_R137_T29TQJ_20250907T150514.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250907T110631_N0511_R137_T30TUP_20250907T150514.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250905T112109_N0511_R037_T30TUP_20250905T133452.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250905T112109_N0511_R037_T29TQJ_20250905T133452.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250902T110619_N0511_R137_T30TUP_20250902T135238.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250729T110641_N0511_R137_T30TUP_20250729T152017.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250727T112119_N0511_R037_T30TUP_20250727T133352.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250729T110641_N0511_R137_T29TQJ_20250729T152017.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250727T112119_N0511_R037_T29TQJ_20250727T133352.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250528T112119_N0511_R037_T30TUP_20250528T133511.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250528T112119_N0511_R037_T29TQJ_20250528T133511.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250525T112131_N0511_R037_T30TUP_20250525T173322.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250525T112131_N0511_R037_T29TQJ_20250525T173322.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250525T110619_N0511_R137_T30TUP_20250525T133604.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250525T110619_N0511_R137_T29TQJ_20250525T133604.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250522T110711_N0511_R137_T30TUP_20250522T150814.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250523T112131_N0511_R037_T30TUP_20250523T134256.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250523T112131_N0511_R037_T29TQJ_20250523T134256.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250518T112119_N0511_R037_T30TUP_20250518T133728.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250518T112119_N0511_R037_T29TQJ_20250518T133728.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250513T112131_N0511_R037_T30TUP_20250513T135636.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250513T112131_N0511_R037_T29TQJ_20250513T135636.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250510T110641_N0511_R137_T29TQJ_20250510T182913.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250508T112119_N0511_R037_T30TUP_20250508T133101.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2A_MSIL2A_20250505T112131_N0511_R037_T30TUP_20250506T030002.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2B_MSIL2A_20250505T110619_N0511_R137_T29TQJ_20250505T134959.zarr β€” Successfully parsed and clipped
βœ… INCLUDE S2C_MSIL2A_20250503T112131_N0511_R037_T30TUP_20250503T165214.zarr β€” Successfully parsed and clipped

βœ“ Included: 109 | βœ— Skipped: 19

Intertidal zone calculation

Water Frequencies

The Scene Classification Layer (SCL)

Every Sentinel-2 L2A product includes a Scene Classification Layer (SCL) β€” a 20-m resolution per-pixel quality and land-cover flag generated by the Sen2Cor atmospheric correction processor. It assigns each pixel one of 12 classes:

SCL value Meaning
0 No data
1 Saturated / defective
2 Dark area pixels
3 Cloud shadows
4 Vegetation
5 Not-vegetated (bare soil/sand)
6 Water
7 Unclassified
8 Cloud (medium probability)
9 Cloud (high probability)
10 Thin cirrus
11 Snow / ice

We focus on classes 4, 5, and 6 (valid clear-sky observations) and use class 6 as the water signal. Other classes indicate cloud, shadow, or sensor issues and are excluded.

Why water frequency reveals the intertidal zone

At any given satellite overpass the estuary is in a particular tidal state: fully exposed, fully submerged, or partially inundated. No single image shows the complete intertidal boundary. However, over many acquisitions that span the full tidal cycle, each pixel accumulates enough observations to estimate the fraction of clear-sky time it spent under water:

\[\text{water\_frequency}(x, y) = \frac{\sum_t \mathbf{1}[\text{SCL}_{t}(x,y) = 6]}{\sum_t \mathbf{1}[\text{SCL}_{t}(x,y) \in \{4, 5, 6\}]}\]

  • Values β‰ˆ 0 β†’ permanently dry land (above high-tide line)
  • Values β‰ˆ 1 β†’ permanently submerged (below low-tide line)
  • Intermediate values (0.1 – 0.95) β†’ intermittently inundated pixels, i.e., the intertidal zone

This approach requires no tidal model, no DEM, and no field calibration: the tidal signal is encoded in the archive itself.

Grid alignment: Because different tiles may be projected in different UTM zones (29N or 30N), we reproject every SCL array onto a common reference grid (datasets[0]) using nearest-neighbour resampling before accumulating counts. This ensures the pixel-wise sums are geographically consistent.

# -----------------------------
# PARAMETERS
# -----------------------------
# Pick a reference dataset for alignment
ref_ds = datasets[0]['scl']

# Reference CRS, transform, and shape
target_crs = ref_ds.rio.crs
ny, nx = ref_ds.sizes['y'], ref_ds.sizes['x']

# Initialize accumulators
water_count = np.zeros((ny, nx), dtype=int)
obs_count = np.zeros((ny, nx), dtype=int)

# -----------------------------
# LOOP OVER DATASETS
# -----------------------------
for ds in datasets:
    if 'scl' not in ds:
        continue

    scl = ds['scl']

    # Flip Y if descending
    if scl.y[0] > scl.y[-1]:
        scl = scl.isel(y=slice(None, None, -1))

    # Ensure CRS is set
    if scl.rio.crs is None:
        scl = scl.rio.write_crs(target_crs, inplace=True)

    # Reproject to reference dataset grid
    scl_rp = scl.rio.reproject_match(
        ref_ds,
        resampling=rasterio.enums.Resampling.nearest
    )

    # Convert to numpy
    scl_vals = scl_rp.values

    # Valid pixels and water mask
    valid = np.isin(scl_vals, [4, 5, 6])  # only clear land or water
    water = (scl_vals == 6)

    # Accumulate counts
    obs_count[valid] += 1
    water_count[valid] += water[valid]

# -----------------------------
# COMPUTE WATER FREQUENCY
# -----------------------------
freq = np.full_like(water_count, fill_value=np.nan, dtype=float)
mask = obs_count > 0
freq[mask] = water_count[mask] / obs_count[mask]

# -----------------------------
# SAVE AS GEO-TIFF
# -----------------------------
out_path = "outputs/water_frequency.tif"

# with rasterio.open(
#     out_path,
#     'w',
#     driver='GTiff',
#     height=ny,
#     width=nx,
#     count=1,
#     dtype='float32',
#     crs=target_crs,
#     transform=ref_ds.rio.transform(),
#     compress='LZW'
# ) as dst:
#     dst.write(freq.astype('float32'), 1)

print(f"Water frequency map saved: {out_path}")

# -----------------------------
# PLOT ON CONTEXTILY BASEMAP
# -----------------------------
# Reproject freq raster to Web Mercator for contextily
_t, _w, _h = calculate_default_transform(
    target_crs, crs_web_mercator, nx, ny, *ref_ds.rio.bounds()
)
freq_wm = np.empty((_h, _w), dtype=np.float32)
rio_reproject(
    source=freq.astype(np.float32),
    destination=freq_wm,
    src_transform=ref_ds.rio.transform(),
    src_crs=target_crs,
    dst_transform=_t,
    dst_crs=crs_web_mercator,
    resampling=Resampling.bilinear,
)
# imshow extent: [left, right, bottom, top]
_ext = [_t.c, _t.c + _t.a * _w, _t.f + _t.e * _h, _t.f]

fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(
    freq_wm, cmap="Blues", vmin=0, vmax=1,
    extent=_ext, origin="upper", alpha=0.7, zorder=2
)
ctx.add_basemap(ax, crs=crs_web_mercator, source=ctx.providers.CartoDB.Positron, zoom=13, zorder=1)
plt.colorbar(im, ax=ax, label="Water frequency", shrink=0.6)
ax.set_title("Water Frequency", fontsize=12, fontweight="bold")
ax.set_axis_off()
plt.tight_layout()
plt.show()
Water frequency map saved: outputs/water_frequency.tif

Intertidal Mask

Choosing the thresholds

We convert the continuous water-frequency raster into a binary intertidal mask using:

\[\text{intertidal}(x, y) = \begin{cases} 1 & \text{if } 0.10 \leq \text{water\_frequency}(x, y) \leq 0.95 \\ 0 & \text{otherwise} \end{cases}\]

The choice of 0.10 and 0.95 is deliberate:

  • Lower bound 0.10: A pixel seen as water in β‰₯10 % of observations has been genuinely inundated β€” not just a noise artefact or a single erroneous SCL classification. Pixels below this threshold are above the high-tide line (supratidal).
  • Upper bound 0.95: A pixel that is water in β‰₯95 % of observations is essentially permanently submerged (subtidal or the main channel). We exclude it because it is not part of the vegetated intertidal habitat we are monitoring.

Adjusting these thresholds lets you expand or contract the intertidal zone definition to match a specific habitat type (e.g. use 0.05–0.50 to focus on the higher, drier saltmarsh fringe).

From raster to vector

The binary mask is vectorised with rasterio.features.shapes, which converts runs of contiguous pixels into GeoJSON polygons. Dissolving all polygons into a single multi-part geometry with geopandas removes interior boundaries and produces a clean habitat layer ready for: - Field survey planning (targeting specific sub-zones) - Overlay with other ecological datasets - Import into GIS software

We threshold the water frequency raster to isolate pixels that are intermittently inundated β€” the classic spectral signature of intertidal habitat. The resulting binary raster is then vectorised using rasterio.features.shapes and dissolved into a single multi-part polygon with geopandas, making it ready for downstream habitat classification or field campaign planning.

Outputs generated by this step

File Format Contents
outputs/intertidal_mask.tif GeoTIFF (uint8) Binary raster: 1 = intertidal, 0 = outside
outputs/intertidal_mask.gpkg GeoPackage (vector) Dissolved polygon of the intertidal zone
(plot) matplotlib figure Mask overlaid on CartoDB basemap 
#-----------------------------
#LOAD WATER FREQUENCY RASTER
#Reproject to ref_ds grid so intertidal_mask always matches scl_arr shape
#-----------------------------
_ny, _nx = ref_ds.sizes["y"], ref_ds.sizes["x"]
_dst_transform = ref_ds.rio.transform()
_dst_crs = ref_ds.rio.crs

with rasterio.open("outputs/water_frequency.tif") as src:
    freq = np.empty((_ny, _nx), dtype=np.float32)

    rio_reproject(
        source=rasterio.band(src, 1),
        destination=freq,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=_dst_transform,
        dst_crs=_dst_crs,
        resampling=Resampling.bilinear,
    )

    transform = _dst_transform
    crs = _dst_crs
    meta = src.meta.copy()

meta.update(height=_ny, width=_nx, transform=transform, crs=crs)

#--------------------

# -----------------------------
# CREATE INTERTIDAL MASK
# -----------------------------
intertidal_mask = np.zeros_like(freq, dtype=np.uint8)
intertidal_mask[(freq >= 0.1) & (freq <= 0.95)] = 1

# -----------------------------
# SAVE MASK AS GEOTIFF
# -----------------------------
# meta.update(dtype=rasterio.uint8, count=1, compress="LZW")

# with rasterio.open("outputs/intertidal_mask.tif", "w", **meta) as dst:
#     dst.write(intertidal_mask, 1)

# print("Intertidal mask raster saved")

# -----------------------------
# RASTER β†’ VECTOR POLYGONS
# -----------------------------
geoms = []

for geom, value in shapes(intertidal_mask, transform=transform):
    if value == 1:
        geoms.append(shape(geom))

gdf = gpd.GeoDataFrame(geometry=geoms, crs=crs)

# -----------------------------
# DISSOLVE INTO SINGLE MULTIPART POLYGON
# -----------------------------
gdf = gdf.dissolve()
gdf = gdf.reset_index(drop=True)

# -----------------------------
# SAVE VECTOR FILE
# -----------------------------
# gdf.to_file("outputs/intertidal_mask.gpkg", driver="GPKG")

# print("Intertidal vector saved: outputs/intertidal_mask.gpkg")

# -----------------------------
# PLOT ON CONTEXTILY BASEMAP
# -----------------------------
gdf_wm = gdf.to_crs(crs_web_mercator)

fig, ax = plt.subplots(figsize=(10, 8))
gdf_wm.plot(ax=ax, color="#1a7abf", alpha=0.6, zorder=2)
ctx.add_basemap(ax, crs=crs_web_mercator, source=ctx.providers.CartoDB.Positron, zoom=13, zorder=1)
ax.set_title("Intertidal Mask", fontsize=12, fontweight="bold")
ax.set_axis_off()
plt.tight_layout()
plt.show()

πŸ’ͺ Now it is your turn

Task 1: Explore Your Own Area of Interest

  • Update bbox and site in the Study Area cell; all downstream cells are bbox-driven. ### Task 2: Vegetation Index
  • Replace the B04/B08 NDVI formula with EVI, SAVI, or a water-quality index such as CHL-Green. ### Task 3: Temporal Analysis Update time_range in the STAC Search cell to restrict or extend the query.

Conclusion

In this notebook we built a fully cloud-native intertidal vegetation monitoring workflow for the Villaviciosa Ramsar estuary using Sentinel-2 L2A data served through the EOPF Zarr Sample Service. The key technical achievements were:

  1. Zero-download analysis β€” the entire workflow ran against cloud object storage via the eopf-zarr xarray engine. No granule was ever saved to disk before processing.
  2. Scalable data discovery β€” pystac_client queried the EOPF STAC catalogue to identify relevant tiles without inspecting any imagery upfront.
  3. Multi-temporal tidal analysis β€” stacking the Scene Classification Layer across ~dozens of acquisitions automatically delineated the intertidal zone from first principles, with no tidal model or DEM required.
  4. Tide-decoupled vegetation monitoring β€” by normalising NDVI-positive pixels against exposed (not total) intertidal area, we removed the tidal confounding effect that would otherwise dominate the time series.
  5. Reproducible and updateable β€” because the EOPF rolling archive is always growing, re-running the notebook at any future date will automatically incorporate new Sentinel-2 acquisitions and refine both the intertidal mask and the NDVI trend.

Ecological findings (preliminary)

The following interpretations are preliminary and will strengthen as the archive grows beyond ~1 year.

  • The water-frequency map clearly distinguishes the main estuarine channel (permanently submerged) from the flanking mudflats and saltmarsh (intermittently inundated).
  • NDVI values above 0.2 on exposed intertidal pixels indicate the presence of actively photosynthesising vegetation (Spartina, macroalgae, seagrass or others).
  • Seasonal differences in vegetated fraction, where sufficient scenes exist, are consistent with the expected temperate phenology (late spring peak, winter minimum).

Limitations & extensions

Limitation Potential improvement
Rolling ~1-year archive Access full historical archive via CDSE or AWS Open Data to extend the time series
Single-site Loop over multiple Ramsar estuaries for a regional monitoring product
SCL-based water detection Replace with a spectral water index (MNDWI, AWEI) computed directly from bands for improved accuracy in turbid estuaries
NDVI cannot distinguish vegetation types NDVI responds to any photosynthesising surface β€” it cannot separate seagrass from macroalgae or other plants. Seagrass is arguably the most ecologically and conservation-critical species in the Villaviciosa intertidal (blue-carbon sink, fish nursery habitat), yet it is spectrally similar to green macroalgae at 10 m. A supervised d spectral classification model trained on multi-band reflectance would enable finer-level mapping and allow seagrass cover to be tracked independently of other vegetation types.
No uncertainty quantification Incorporate field data to derive confidence intervals

What’s next?

In the following notebook we will dive into Sentinel-2 L2A application for cloud probability assessment based on the cld layer.

This notebook was developed as a contribution to the EOPF Zarr Community Notebook Competition.