Search NASA OPERA granules¶

# %pip install geolibre geopandas requests matplotlib

import re

import geopandas as gpd
import pandas as pd
import requests
from geolibre import Map
from shapely.geometry import Polygon, box

OPERA_PRODUCTS = {
    "OPERA_L3_DSWX-HLS_V1": "Dynamic Surface Water Extent from HLS",
    "OPERA_L3_DSWX-S1_V1": "Dynamic Surface Water Extent from Sentinel-1",
    "OPERA_L3_DIST-ALERT-HLS_V1": "Land Surface Disturbance Alert from HLS",
    "OPERA_L3_DIST-ANN-HLS_V1": "Land Surface Disturbance Annual from HLS",
    "OPERA_L2_RTC-S1_V1": "Radiometric Terrain Corrected SAR Backscatter from Sentinel-1",
    "OPERA_L2_RTC-S1-STATIC_V1": "RTC-S1 Static Layers",
    "OPERA_L2_CSLC-S1_V1": "Coregistered Single-Look Complex from Sentinel-1",
    "OPERA_L2_CSLC-S1-STATIC_V1": "CSLC-S1 Static Layers",
}
pd.Series(OPERA_PRODUCTS, name="title").rename_axis("short_name").to_frame()

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)

bbox = [-0.55, 39.2, -0.25, 39.5]  # west, south, east, north
start_date, end_date = "2024-10-25", "2024-11-10"

dswx_hls = search_granules("OPERA_L3_DSWX-HLS_V1", bbox, start_date, end_date)
print(f"Found {len(dswx_hls)} granules")
dswx_hls.drop(columns=["urls", "geometry"]).head(10)

dswx_hls.loc[0, "bands"].split(",")

dswx_hls.loc[0, "urls"]["B01_WTR"]

m = Map(center=(-0.4, 39.35), zoom=8, height="600px")
m

m.add_gdf(
    dswx_hls.drop(columns="urls"),
    name="DSWx-HLS footprints",
    fillColor="#3b82f6",
    fillOpacity=0.1,
    strokeColor="#1d4ed8",
    strokeWidth=1,
)
search_area = m.add_geojson(
    box(*bbox).__geo_interface__,
    name="Search area",
    fillOpacity=0,
    strokeColor="#ef4444",
    strokeWidth=2,
)

products = {
    "OPERA_L3_DSWX-S1_V1": "#06b6d4",
    "OPERA_L2_RTC-S1_V1": "#a855f7",
    "OPERA_L3_DIST-ALERT-HLS_V1": "#f97316",
}

results = {"OPERA_L3_DSWX-HLS_V1": dswx_hls}
for short_name, color in products.items():
    gdf = search_granules(short_name, bbox, start_date, end_date)
    results[short_name] = gdf
    if len(gdf):
        m.add_gdf(
            gdf.drop(columns="urls"),
            name=f"{short_name} footprints",
            fillColor=color,
            fillOpacity=0.05,
            strokeColor=color,
            strokeWidth=1,
        )

summary = pd.DataFrame(
    {
        "granules": {k: len(v) for k, v in results.items()},
        "first": {k: v["begin_date"].min()[:10] for k, v in results.items() if len(v)},
        "last": {k: v["begin_date"].max()[:10] for k, v in results.items() if len(v)},
    }
)
summary

try:
    drawn = m.get_drawn_features(as_gdf=True)
except TimeoutError:
    drawn = []  # the map is not displayed, e.g. when the notebook runs headless

if len(drawn):
    drawn_bbox = [round(v, 4) for v in drawn.total_bounds]
    drawn_hits = search_granules(
        "OPERA_L3_DSWX-HLS_V1", drawn_bbox, start_date, end_date
    )
    print(f"{len(drawn_hits)} DSWx-HLS granules in {drawn_bbox}")
else:
    print("Draw a rectangle on the map first.")

dswx_hls.drop(columns="urls").to_file("dswx_hls_footprints.geojson", driver="GeoJSON")