Map a flood with OPERA DSWx-S1¶

# %pip install geolibre geopandas requests matplotlib

import json
import re

import geopandas as gpd
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import pandas as pd
import requests
from geolibre import Map
from shapely.geometry import Polygon, box

CMR_GRANULES_URL = "https://cmr.earthdata.nasa.gov/search/granules.umm_json"
BAND_PATTERN = re.compile(r"_(B\d{2}_[A-Za-z0-9-]+|[VH]{2}|[A-Z][A-Za-z0-9-]+)\.tiff?$")


def granule_geometry(umm):
    """Build a footprint polygon from a UMM-G granule record.

    Args:
        umm: The ``umm`` block of a CMR UMM-G granule item.

    Returns:
        A shapely Polygon, or None when the record has no usable geometry.
    """
    geometry = umm.get("SpatialExtent", {}).get("HorizontalSpatialDomain", {})
    geometry = geometry.get("Geometry", {})
    rects = geometry.get("BoundingRectangles") or []
    polygons = geometry.get("GPolygons") or []
    if rects:
        r = rects[0]
        return box(
            r["WestBoundingCoordinate"],
            r["SouthBoundingCoordinate"],
            r["EastBoundingCoordinate"],
            r["NorthBoundingCoordinate"],
        )
    if polygons:
        points = polygons[0]["Boundary"]["Points"]
        return Polygon([(p["Longitude"], p["Latitude"]) for p in points])
    return None


def search_granules(short_name, bbox, start_date, end_date, count=100):
    """Search NASA CMR for OPERA granules (no Earthdata login required).

    Args:
        short_name: OPERA collection short name, e.g. ``"OPERA_L3_DSWX-HLS_V1"``.
        bbox: Search area as ``[west, south, east, north]`` in degrees.
        start_date: First day of the search window, ``"YYYY-MM-DD"``.
        end_date: Last day of the search window, ``"YYYY-MM-DD"``.
        count: Maximum number of granules to return.

    Returns:
        A GeoDataFrame with one footprint per granule, sorted by date.
    """
    params = {
        "short_name": short_name,
        "bounding_box": ",".join(str(v) for v in bbox),
        "temporal": f"{start_date}T00:00:00Z,{end_date}T23:59:59Z",
        "page_size": count,
    }
    response = requests.get(CMR_GRANULES_URL, params=params, timeout=60)
    response.raise_for_status()

    rows = []
    for item in response.json().get("items", []):
        umm = item["umm"]
        dates = umm.get("TemporalExtent", {}).get("RangeDateTime", {})
        # CMR lists each file twice (HTTPS and S3); keep the HTTPS copy.
        urls = {}
        for link in umm.get("RelatedUrls", []):
            url = link.get("URL", "")
            match = BAND_PATTERN.search(url)
            if url.startswith("https://") and match:
                urls.setdefault(match.group(1), url)
        rows.append(
            {
                "granule_ur": umm["GranuleUR"],
                "concept_id": item["meta"]["collection-concept-id"],
                "begin_date": dates.get("BeginningDateTime", ""),
                "bands": ",".join(urls),
                "urls": urls,
                "geometry": granule_geometry(umm),
            }
        )

    columns = ["granule_ur", "concept_id", "begin_date", "bands", "urls", "geometry"]
    gdf = gpd.GeoDataFrame(rows, columns=columns, geometry="geometry", crs="EPSG:4326")
    return gdf.sort_values("begin_date").reset_index(drop=True)

TITILER_CMR = "https://titiler-cmr.opengeos.org"


def titiler_params(granule, band, assets_regex, **render):
    """Build the titiler-cmr query that pins one granule and one band.

    Args:
        granule: A row of the GeoDataFrame returned by ``search_granules``.
        band: Band (asset) name, e.g. ``"B01_WTR"`` or ``"VV"``.
        assets_regex: Regex titiler-cmr uses to find band names in file names.
        **render: Extra titiler parameters such as ``rescale``,
            ``colormap_name``, ``colormap`` (a dict), or ``expression``.

    Returns:
        A dict of query parameters.
    """
    params = {
        "collection_concept_id": granule["concept_id"],
        "granule_ur": granule["granule_ur"],
        "assets": band,
        "assets_regex": assets_regex,
    }
    for key, value in render.items():
        if value is None:
            continue
        if key == "colormap":
            value = json.dumps(value)
        params[key] = value
    if "expression" in params:
        # Expressions refer to the requested asset as b1.
        params["asset_as_band"] = "true"
    return params


def opera_tiles(granule, band, assets_regex, **render):
    """Ask titiler-cmr for an XYZ tile URL that renders one OPERA granule.

    titiler-cmr reads the Cloud-Optimized GeoTIFF and handles Earthdata
    authentication server-side, so no NASA login is needed here.

    Args:
        granule: A row of the GeoDataFrame returned by ``search_granules``.
        band: Band (asset) name to render.
        assets_regex: Regex titiler-cmr uses to find band names in file names.
        **render: Extra titiler parameters, see ``titiler_params``.

    Returns:
        A ``(tile_url, bounds)`` tuple for ``Map.add_tile_layer``.
    """
    url = f"{TITILER_CMR}/rasterio/WebMercatorQuad/tilejson.json"
    params = titiler_params(granule, band, assets_regex, **render)
    response = requests.get(url, params=params, timeout=120)
    response.raise_for_status()
    tilejson = response.json()
    return tilejson["tiles"][0], tilejson.get("bounds")

def class_areas(granule, aoi, band="B01_WTR", assets_regex=r"B[0-9]{2}_[A-Z]+"):
    """Count the area of each class of a categorical OPERA band inside an AOI.

    OPERA rasters use a 30 m grid, and titiler-cmr computes statistics in the
    granule's native projection, so each pixel is 900 square meters.

    Args:
        granule: A row of the GeoDataFrame returned by ``search_granules``.
        aoi: A shapely geometry in EPSG:4326.
        band: Band (asset) name.
        assets_regex: Regex titiler-cmr uses to find band names in file names.

    Returns:
        A dict mapping class value to area in square kilometers.
    """
    params = titiler_params(granule, band, assets_regex, categorical="true")
    feature = {"type": "Feature", "properties": {}, "geometry": aoi.__geo_interface__}
    response = requests.post(
        f"{TITILER_CMR}/rasterio/statistics", params=params, json=feature, timeout=180
    )
    response.raise_for_status()
    counts, classes = response.json()["properties"]["statistics"]["b1"]["histogram"]
    return {int(c): n * 900 / 1e6 for c, n in zip(classes, counts)}

aoi_bbox = [-0.45, 39.25, -0.30, 39.45]
aoi = box(*aoi_bbox)

m = Map(center=(-0.375, 39.35), zoom=10, height="650px")
m.add_geojson(
    aoi.__geo_interface__,
    name="Area of interest",
    fillOpacity=0,
    strokeColor="#ef4444",
    strokeWidth=2,
)
m

scenes = search_granules("OPERA_L3_DSWX-S1_V1", aoi_bbox, "2024-10-20", "2024-11-10")
scenes = scenes[scenes["granule_ur"].str.contains("_T30SYJ_")].copy()
scenes["date"] = pd.to_datetime(scenes["begin_date"]).dt.date
scenes = scenes.drop_duplicates("date").reset_index(drop=True)
scenes[["date", "granule_ur"]]

DSWX_S1_COLORMAP = {
    "0": [0, 0, 0, 0],  # not water
    "1": [0, 0, 255, 255],  # open water
    "3": [0, 168, 132, 255],  # inundated vegetation
    "250": [0, 0, 0, 0],  # HAND masked
    "251": [0, 0, 0, 0],  # layover/shadow masked
    "255": [0, 0, 0, 0],  # fill
}


def add_scene(granule, name):
    """Add one DSWx-S1 scene to the map as a water-classification layer.

    Args:
        granule: A row of the GeoDataFrame returned by ``search_granules``.
        name: Layer name shown in the layer panel.

    Returns:
        The id of the added layer.
    """
    url, bounds = opera_tiles(
        granule, "B01_WTR", r"B[0-9]{2}_[A-Z]+", colormap=DSWX_S1_COLORMAP
    )
    return m.add_tile_layer(url, name=name, bounds=bounds, attribution="NASA OPERA")


dates = scenes["date"].astype(str)
before = scenes[dates == "2024-10-25"].iloc[0]
after = scenes[dates == "2024-10-31"].iloc[0]

before_layer = add_scene(before, "DSWx-S1 before (2024-10-25)")
after_layer = add_scene(after, "DSWx-S1 after (2024-10-31)")
m.split_map(left_layers=before_layer, right_layers=after_layer)
m.add_legend(
    "DSWx-S1",
    legend_dict={"Open water": "#0000ff", "Inundated vegetation": "#00a884"},
)

CLASS_NAMES = {1: "Open water", 3: "Inundated vegetation"}

before_areas = class_areas(before, aoi)
after_areas = class_areas(after, aoi)

comparison = pd.DataFrame(
    {
        "before (km²)": {CLASS_NAMES[c]: before_areas.get(c, 0) for c in CLASS_NAMES},
        "after (km²)": {CLASS_NAMES[c]: after_areas.get(c, 0) for c in CLASS_NAMES},
    }
)
comparison["change (km²)"] = comparison["after (km²)"] - comparison["before (km²)"]
comparison.round(1)

records = []
for _, scene in scenes.iterrows():
    areas = class_areas(scene, aoi)
    records.append(
        {
            "date": scene["date"],
            "Open water": areas.get(1, 0),
            "Inundated vegetation": areas.get(3, 0),
        }
    )
timeseries = pd.DataFrame(records)
timeseries["date"] = pd.to_datetime(timeseries["date"])
timeseries = timeseries.set_index("date")
timeseries.round(1)

fig, ax = plt.subplots(figsize=(9, 4))
timeseries.plot(ax=ax, marker="o", color=["#1d4ed8", "#0d9488"])
ax.axvline(
    pd.Timestamp("2024-10-29"), color="#dc2626", linestyle="--", label="DANA flood"
)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %d"))
ax.set_ylabel("Area (km²)")
ax.set_xlabel("")
ax.set_title("DSWx-S1 surface water in the AOI")
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()