From 56205877e7a831276d0e2fbbb8ead7bc96f2b7aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Tue, 26 May 2026 10:56:02 +0100 Subject: [PATCH 01/19] feat(stac): add build_s1_rtc_stac_item and generate-stac-s1 CLI command Adds eopf_geozarr.stac.s1_rtc.build_s1_rtc_stac_item, which reads a consolidated S1 GRD RTC Zarr V3 store and returns a pystac.Item with SAR/SAT/projection extensions, WGS84 bbox derived via pyproj, and vv/vh asset sub-paths. Ascending orbit is preferred when both are present. Also adds the generate-stac-s1 CLI subcommand that prints the item as JSON. 8 unit tests cover all acceptance criteria from the plan. Co-Authored-By: Claude Sonnet 4.6 --- pyproject.toml | 1 + src/eopf_geozarr/cli.py | 33 ++++++ src/eopf_geozarr/stac/__init__.py | 1 + src/eopf_geozarr/stac/s1_rtc.py | 147 +++++++++++++++++++++++ tests/test_s1_stac.py | 191 ++++++++++++++++++++++++++++++ 5 files changed, 373 insertions(+) create mode 100644 src/eopf_geozarr/stac/__init__.py create mode 100644 src/eopf_geozarr/stac/s1_rtc.py create mode 100644 tests/test_s1_stac.py diff --git a/pyproject.toml b/pyproject.toml index 47a8017e..c55fb102 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,6 +41,7 @@ dependencies = [ "s3fs>=2024.6.0", "boto3>=1.34.0", "pyproj>=3.7.0", + "pystac>=1.8.0", "structlog>=25.5.0", ] diff --git a/src/eopf_geozarr/cli.py b/src/eopf_geozarr/cli.py index ea4bc867..7735cac5 100755 --- a/src/eopf_geozarr/cli.py +++ b/src/eopf_geozarr/cli.py @@ -1108,6 +1108,20 @@ def consolidate_s1_command(args: argparse.Namespace) -> None: sys.exit(1) +def generate_stac_s1_command(args: argparse.Namespace) -> None: + """Build and print a STAC item for an S1 GRD RTC Zarr store.""" + import json + + from .stac.s1_rtc import build_s1_rtc_stac_item + + try: + item = build_s1_rtc_stac_item(args.store, args.collection) + print(json.dumps(item.to_dict(), indent=2)) + except Exception as e: + log.exception("❌ Error generating STAC item", error=str(e)) + sys.exit(1) + + def create_parser() -> argparse.ArgumentParser: """ Create the argument parser for the CLI. @@ -1238,6 +1252,7 @@ def create_parser() -> argparse.ArgumentParser: # Add S1 ingestion commands add_s1_ingestion_commands(subparsers) + add_s1_stac_commands(subparsers) return parser @@ -1304,6 +1319,24 @@ def add_s1_ingestion_commands(subparsers: argparse._SubParsersAction) -> None: cons_parser.set_defaults(func=consolidate_s1_command) +def add_s1_stac_commands(subparsers: argparse._SubParsersAction) -> None: + """Add S1 GRD RTC STAC builder command to CLI parser.""" + stac_parser = subparsers.add_parser( + "generate-stac-s1", + help="Build and print a STAC item for an S1 GRD RTC Zarr store", + ) + stac_parser.add_argument( + "--store", type=str, required=True, help="Path or s3:// URI to the Zarr store" + ) + stac_parser.add_argument( + "--collection", + type=str, + default="sentinel-1-grd-rtc-staging", + help="STAC collection ID (default: sentinel-1-grd-rtc-staging)", + ) + stac_parser.set_defaults(func=generate_stac_s1_command) + + def add_s2_optimization_commands(subparsers: argparse._SubParsersAction) -> None: """Add S2 optimization commands to CLI parser.""" diff --git a/src/eopf_geozarr/stac/__init__.py b/src/eopf_geozarr/stac/__init__.py new file mode 100644 index 00000000..3ce158bd --- /dev/null +++ b/src/eopf_geozarr/stac/__init__.py @@ -0,0 +1 @@ +"""STAC item builders for eopf-geozarr Zarr stores.""" diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py new file mode 100644 index 00000000..3c57c567 --- /dev/null +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -0,0 +1,147 @@ +"""STAC item builder for S1 GRD RTC Zarr V3 stores.""" + +from __future__ import annotations + +import datetime as dt +from pathlib import Path +from typing import cast + +import numpy as np +import pyproj +import pystac +import zarr + +SAR_EXT = "https://stac-extensions.github.io/sar/v1.0.0/schema.json" +SAT_EXT = "https://stac-extensions.github.io/sat/v1.0.0/schema.json" +PROJ_EXT = "https://stac-extensions.github.io/projection/v2.0.0/schema.json" + +_ORBIT_PREFERENCE = ("ascending", "descending") + + +def _utm_to_wgs84( + proj_code: str, utm_bbox: list[float] +) -> tuple[float, float, float, float]: + """Convert UTM [xmin, ymin, xmax, ymax] to WGS84 (west, south, east, north).""" + xmin, ymin, xmax, ymax = utm_bbox + transformer = pyproj.Transformer.from_crs(proj_code, "EPSG:4326", always_xy=True) + xs = [xmin, xmax, xmin, xmax] + ys = [ymin, ymin, ymax, ymax] + lons, lats = transformer.transform(xs, ys) + return min(lons), min(lats), max(lons), max(lats) + + +def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: + """Build a STAC item from a consolidated S1 GRD RTC Zarr V3 store. + + Parameters + ---------- + zarr_store: + Local path or ``s3://`` URI to the Zarr store. + collection_id: + STAC collection ID to attach to the item. + + Returns + ------- + pystac.Item + + Raises + ------ + ValueError + If the store contains no acquisitions. + """ + tile_id = Path(zarr_store).name.removeprefix("s1-grd-rtc-").removesuffix(".zarr") + + root = zarr.open_consolidated(zarr_store, zarr_format=3) + + all_times_ns: list[int] = [] + wgs84_bboxes: list[tuple[float, float, float, float]] = [] + preferred_orbit: str | None = None + + for orbit_dir in _ORBIT_PREFERENCE: + if orbit_dir not in root: + continue + og = cast(zarr.Group, root[orbit_dir]) + attrs = dict(og.attrs) + proj_code = cast(str, attrs["proj:code"]) + utm_bbox = cast(list[float], attrs["spatial:bbox"]) + + r10m = cast(zarr.Group, og["r10m"]) + times = np.array(cast(zarr.Array, r10m["time"])).tolist() + if not times: + continue + + all_times_ns.extend(times) + wgs84_bboxes.append(_utm_to_wgs84(proj_code, utm_bbox)) + if preferred_orbit is None: + preferred_orbit = orbit_dir + + if not all_times_ns: + raise ValueError(f"No acquisitions found in Zarr store: {zarr_store}") + + # Temporal range + start_dt = dt.datetime.fromtimestamp(min(all_times_ns) / 1e9, tz=dt.UTC) + end_dt = dt.datetime.fromtimestamp(max(all_times_ns) / 1e9, tz=dt.UTC) + + # WGS84 bbox union across all present orbit directions + west = min(b[0] for b in wgs84_bboxes) + south = min(b[1] for b in wgs84_bboxes) + east = max(b[2] for b in wgs84_bboxes) + north = max(b[3] for b in wgs84_bboxes) + wgs84_bbox = [west, south, east, north] + + geometry = { + "type": "Polygon", + "coordinates": [ + [[west, south], [east, south], [east, north], [west, north], [west, south]] + ], + } + + # preferred_orbit is guaranteed non-None here (ValueError raised above if no acquisitions) + assert preferred_orbit is not None + preferred_proj_code = cast(str, dict(cast(zarr.Group, root[preferred_orbit]).attrs)["proj:code"]) + + item = pystac.Item( + id=f"s1-rtc-{tile_id}", + geometry=geometry, + bbox=wgs84_bbox, + datetime=None, + properties={ + "start_datetime": start_dt.isoformat(), + "end_datetime": end_dt.isoformat(), + # SAR extension + "sar:instrument_mode": "IW", + "sar:frequency_band": "C", + "sar:center_frequency": 5.405, + "sar:polarizations": ["VV", "VH"], + "sar:product_type": "GRD", + # SAT extension + "sat:orbit_state": preferred_orbit, + # Projection extension + "proj:code": preferred_proj_code, + }, + stac_extensions=[SAR_EXT, SAT_EXT, PROJ_EXT], + collection=collection_id, + ) + + store_str = str(zarr_store) + item.add_asset( + "zarr-store", + pystac.Asset( + href=store_str, + media_type="application/vnd.zarr; version=3", + roles=["data"], + title="S1 GRD RTC Zarr store", + ), + ) + for pol in ("vv", "vh"): + item.add_asset( + pol, + pystac.Asset( + href=f"{store_str}/{preferred_orbit}/r10m/{pol}", + media_type="application/vnd.zarr; version=3", + roles=["data"], + title=f"{pol.upper()} polarisation (native 10 m)", + ), + ) + + return item diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py new file mode 100644 index 00000000..a67a9c87 --- /dev/null +++ b/tests/test_s1_stac.py @@ -0,0 +1,191 @@ +"""Tests for build_s1_rtc_stac_item — STAC item builder for S1 GRD RTC Zarr stores.""" + +from __future__ import annotations + +import datetime as dt +from typing import TYPE_CHECKING + +import numpy as np + +if TYPE_CHECKING: + from pathlib import Path +import pytest +import zarr + +from eopf_geozarr.stac.s1_rtc import build_s1_rtc_stac_item + +# ============================================================================= +# Constants +# ============================================================================= + +CRS = "EPSG:32631" +UTM_BBOX = [300000.0, 4900000.0, 400000.0, 5000000.0] # [xmin, ymin, xmax, ymax] + +# Nanoseconds since epoch for two acquisitions +T1_NS = int(np.datetime64("2023-01-15T06:12:34", "ns").astype(np.int64)) +T2_NS = int(np.datetime64("2023-01-27T06:12:35", "ns").astype(np.int64)) + + +# ============================================================================= +# Fixture helper +# ============================================================================= + + +def _make_s1_store( + tmp_path: Path, + orbits: dict[str, list[tuple[int, str]]], + tile_id: str = "31TCH", + crs: str = CRS, + utm_bbox: list[float] | None = None, +) -> Path: + """Create a minimal consolidated S1 Zarr store. + + ``orbits`` maps orbit_direction -> list of (time_ns, platform) tuples. + Creates only the attrs and coordinate arrays needed by build_s1_rtc_stac_item. + """ + if utm_bbox is None: + utm_bbox = UTM_BBOX + store_path = tmp_path / f"s1-grd-rtc-{tile_id}.zarr" + root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + for orbit_dir, acquisitions in orbits.items(): + og = root.create_group(orbit_dir) + og.attrs.update({"proj:code": crs, "spatial:bbox": utm_bbox}) + r10m = og.create_group("r10m") + times = np.array([t for t, _ in acquisitions], dtype="int64") + platforms = np.array([p for _, p in acquisitions], dtype=" None: + """Item id must be s1-rtc-{tile_id} derived from the store basename.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + assert item.id == "s1-rtc-31TCH" + + +def test_temporal_range_min_max(tmp_path: Path) -> None: + """start_datetime/end_datetime must span the full time range across all acquisitions.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A"), (T2_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + start = dt.datetime.fromisoformat(item.properties["start_datetime"]) + end = dt.datetime.fromisoformat(item.properties["end_datetime"]) + + expected_start = dt.datetime(2023, 1, 15, 6, 12, 34, tzinfo=dt.UTC) + expected_end = dt.datetime(2023, 1, 27, 6, 12, 35, tzinfo=dt.UTC) + + assert abs((start - expected_start).total_seconds()) < 1 + assert abs((end - expected_end).total_seconds()) < 1 + assert item.datetime is None + + +def test_bbox_wgs84_from_utm(tmp_path: Path) -> None: + """UTM bbox must be converted to WGS84 and stored as item bbox.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + west, south, east, north = item.bbox # type: ignore[misc] + # EPSG:32631 [300000,4900000,400000,5000000] -> approx 0.46E-1.75E, 44.2N-45.1N + assert 0.0 < west < 1.0 + assert 44.0 < south < 45.0 + assert 1.0 < east < 2.0 + assert 45.0 < north < 46.0 + + +def test_both_orbits_bbox_union(tmp_path: Path) -> None: + """When ascending and descending are both present, the WGS84 bbox is the union.""" + # Give ascending a different UTM bbox (shifted east) + store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + + for orbit_dir, bbox in [ + ("descending", [300000.0, 4900000.0, 400000.0, 5000000.0]), + ("ascending", [400000.0, 4900000.0, 500000.0, 5000000.0]), + ]: + og = root.create_group(orbit_dir) + og.attrs.update({"proj:code": CRS, "spatial:bbox": bbox}) + r10m = og.create_group("r10m") + t_arr = r10m.create_array("time", shape=(1,), dtype="int64", chunks=(512,)) + t_arr[:] = [T1_NS] + p_arr = r10m.create_array("platform", shape=(1,), dtype=" 2.5 # right edge from ascending (shifted ~1° further east) + + +def test_ascending_preferred_for_assets(tmp_path: Path) -> None: + """When both orbits present, vv/vh assets must point to the ascending group.""" + store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + for orbit_dir in ("descending", "ascending"): + og = root.create_group(orbit_dir) + og.attrs.update({"proj:code": CRS, "spatial:bbox": UTM_BBOX}) + r10m = og.create_group("r10m") + t_arr = r10m.create_array("time", shape=(1,), dtype="int64", chunks=(512,)) + t_arr[:] = [T1_NS] + p_arr = r10m.create_array("platform", shape=(1,), dtype=" None: + """A store with an orbit group but no acquisitions must raise ValueError.""" + store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + og = root.create_group("descending") + og.attrs.update({"proj:code": CRS, "spatial:bbox": UTM_BBOX}) + r10m = og.create_group("r10m") + t_arr = r10m.create_array("time", shape=(0,), dtype="int64", chunks=(512,)) + p_arr = r10m.create_array("platform", shape=(0,), dtype=" None: + """zarr-store href = store URI; vv/vh hrefs = {store}/{orbit}/r10m/{pol}.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + store_str = str(store) + assert item.assets["zarr-store"].href == store_str + assert item.assets["vv"].href == f"{store_str}/descending/r10m/vv" + assert item.assets["vh"].href == f"{store_str}/descending/r10m/vh" + + +def test_sar_extension_fields(tmp_path: Path) -> None: + """SAR extension fields must be set with correct values for S1 IW GRD.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + props = item.properties + assert props["sar:instrument_mode"] == "IW" + assert props["sar:frequency_band"] == "C" + assert props["sar:center_frequency"] == pytest.approx(5.405) + assert props["sar:polarizations"] == ["VV", "VH"] + assert props["sar:product_type"] == "GRD" + + sar_ext_uri = "https://stac-extensions.github.io/sar/v1.0.0/schema.json" + assert sar_ext_uri in item.stac_extensions From 56160cd3b54198ddd35a60732b819f773f2b821f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Fri, 29 May 2026 11:05:49 +0100 Subject: [PATCH 02/19] fix: support s3:// URIs in S1Tiling discovery and ingest (#175) Four functions in s1_ingest.py silently failed for s3:// paths: - discover_s1tiling_acquisitions / discover_s1tiling_conditions used pathlib.Path.glob(), which normalises s3:// to s3:/ and returns 0 matches - ingest_s1tiling_acquisition / ingest_s1tiling_conditions coerced paths via Path(), corrupting s3:// URIs, and called .exists() which always returns False for S3 (existence checks now use path_exists() from fs_utils) Adds three private helpers: - _list_tifs(): uses s3fs.S3FileSystem.glob for s3:// prefixes - _coerce_input_path(): preserves str for s3://, returns Path otherwise - _rasterio_env(): rasterio.Env(AWSSession) context for s3:// paths; strips https:// from AWS_ENDPOINT_URL since GDAL expects hostname only Tested: rasterio.open("s3://...") confirmed working on OVH S3 with AWS_S3_ENDPOINT or AWS_ENDPOINT_URL set; hostname-only required by GDAL. Refs #139 Co-authored-by: Claude Sonnet 4.6 --- src/eopf_geozarr/conversion/s1_ingest.py | 99 ++++++++++++++++++------ tests/test_s1_rtc_ingest.py | 24 ++++++ 2 files changed, 101 insertions(+), 22 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 8109e8a7..938c6f40 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -130,7 +130,7 @@ def extract_geotiff_metadata(path: str | Path) -> S1TilingMetadata: If critical tags (ACQUISITION_DATETIME, ORBIT_NUMBER, RELATIVE_ORBIT_NUMBER, FLYING_UNIT_CODE) are missing. """ - with rasterio.open(str(path)) as src: + with _rasterio_env(path), rasterio.open(str(path)) as src: tags = src.tags() t = src.transform spatial_transform = [t.a, t.b, t.c, t.d, t.e, t.f] @@ -470,13 +470,13 @@ def ingest_s1tiling_acquisition( ValueError If the GeoTIFF CRS or shape does not match the existing store. """ - vv_path = Path(vv_path) - vh_path = Path(vh_path) - border_mask_path = Path(border_mask_path) + vv_path = _coerce_input_path(vv_path) + vh_path = _coerce_input_path(vh_path) + border_mask_path = _coerce_input_path(border_mask_path) store_path = Path(store_path) for p in [vv_path, vh_path, border_mask_path]: - if not p.exists(): + if not _input_path_exists(p): raise FileNotFoundError(f"GeoTIFF not found: {p}") # Extract metadata from VV file @@ -583,12 +583,13 @@ def ingest_s1tiling_acquisition( orbit = root[orbit_direction] # Read GeoTIFF pixel data - with rasterio.open(str(vv_path)) as src: - vv_data = src.read(1) - with rasterio.open(str(vh_path)) as src: - vh_data = src.read(1) - with rasterio.open(str(border_mask_path)) as src: - mask_data = src.read(1).astype(np.uint8) + with _rasterio_env(vv_path): + with rasterio.open(str(vv_path)) as src: + vv_data = src.read(1) + with rasterio.open(str(vh_path)) as src: + vh_data = src.read(1) + with rasterio.open(str(border_mask_path)) as src: + mask_data = src.read(1).astype(np.uint8) log.info( "GeoTIFF read complete", @@ -665,6 +666,62 @@ def consolidate_s1_store(store_path: str | Path, orbit_direction: str) -> None: ) +# ============================================================================= +# S3 / local filesystem helpers +# ============================================================================= + + +def _list_tifs(input_dir: str | Path) -> list[str | Path]: + """List *.tif files; supports local paths and s3:// URIs.""" + s = str(input_dir).rstrip("/") + if s.startswith("s3://"): + import s3fs as _s3 + + fs = _s3.S3FileSystem() + bucket_key = s[len("s3://") :] + return [f"s3://{p}" for p in sorted(fs.glob(f"{bucket_key}/*.tif"))] + return sorted(Path(input_dir).glob("*.tif")) + + +def _coerce_input_path(p: str | Path) -> str | Path: + """Return str for s3:// URIs (preserves double-slash); Path otherwise.""" + s = str(p) + return s if s.startswith("s3://") else Path(s) + + +def _input_path_exists(p: str | Path) -> bool: + """Existence check for both local Path and s3:// URI.""" + s = str(p) + if s.startswith("s3://"): + import s3fs as _s3 + + return _s3.S3FileSystem().exists(s.removeprefix("s3://")) + return Path(p).exists() + + +def _rasterio_env(path: str | Path): # type: ignore[return] + """rasterio.Env context for S3 paths; no-op context manager for local paths. + + rasterio 1.5 passes endpoint_url verbatim as GDAL's AWS_S3_ENDPOINT. + GDAL expects hostname only (no scheme), so strip https:// from + AWS_ENDPOINT_URL if that is what the environment provides. + """ + import contextlib + import os + + if not str(path).startswith("s3://"): + return contextlib.nullcontext() + + import boto3 + import rasterio + from rasterio.session import AWSSession + + raw = os.environ.get("AWS_S3_ENDPOINT") or os.environ.get("AWS_ENDPOINT_URL", "") + endpoint = raw.split("://", 1)[-1] if "://" in raw else raw + session = boto3.Session() + return rasterio.Env(AWSSession(session, endpoint_url=endpoint or None)) + + # ============================================================================= # File Discovery # ============================================================================= @@ -689,12 +746,11 @@ def discover_s1tiling_acquisitions(input_dir: str | Path) -> list[dict]: Logs warnings for incomplete acquisitions (missing polarisation or mask files). """ - input_dir = Path(input_dir) - files = sorted(input_dir.glob("*.tif")) + files = _list_tifs(input_dir) groups: dict[tuple, dict] = {} for f in files: - parsed = parse_s1tiling_filename(f.name) + parsed = parse_s1tiling_filename(Path(str(f)).name) if parsed is None: continue @@ -792,8 +848,8 @@ def ingest_s1tiling_conditions( ("incidence_angle", incidence_angle_path), ]: if path is not None: - p = Path(path) - if not p.exists(): + p = _coerce_input_path(path) + if not _input_path_exists(p): raise FileNotFoundError(f"Condition GeoTIFF not found: {p}") condition_inputs.append((label, p)) @@ -817,7 +873,7 @@ def ingest_s1tiling_conditions( # Read reference metadata from the first condition file ref_label, ref_path = condition_inputs[0] - with rasterio.open(str(ref_path)) as src: + with _rasterio_env(ref_path), rasterio.open(str(ref_path)) as src: ref_crs = str(src.crs) t = src.transform ref_transform = [t.a, t.b, t.c, t.d, t.e, t.f] @@ -849,7 +905,7 @@ def ingest_s1tiling_conditions( for label, cond_path in condition_inputs: array_name = f"{label}_{orbit_suffix}" - with rasterio.open(str(cond_path)) as src: + with _rasterio_env(cond_path), rasterio.open(str(cond_path)) as src: data = src.read(1).astype(np.float32) h, w = data.shape @@ -901,12 +957,11 @@ def discover_s1tiling_conditions(input_dir: str | Path) -> list[dict]: Groups by (tile, orbit). """ - input_dir = Path(input_dir) - files = sorted(input_dir.glob("*.tif")) + files = _list_tifs(input_dir) groups: dict[tuple[str, str], dict] = {} for f in files: - m = S1TILING_GAMMA_AREA_PATTERN.match(f.name) + m = S1TILING_GAMMA_AREA_PATTERN.match(Path(str(f)).name) if m: tile = m.group("tile") orbit = m.group("orbit") @@ -916,7 +971,7 @@ def discover_s1tiling_conditions(input_dir: str | Path) -> list[dict]: groups[key]["gamma_area"] = f continue - m = S1TILING_LIA_PATTERN.match(f.name) + m = S1TILING_LIA_PATTERN.match(Path(str(f)).name) if m: tile = m.group("tile") orbit = m.group("orbit") diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py index c5e4b538..e404d642 100644 --- a/tests/test_s1_rtc_ingest.py +++ b/tests/test_s1_rtc_ingest.py @@ -4,6 +4,7 @@ from math import ceil from pathlib import Path +from unittest.mock import patch import numpy as np import pytest @@ -525,6 +526,20 @@ def test_resolves_masked_multiframe_stamp_from_tag(self, tmp_path: Path) -> None for k in ("vv", "vh", "vv_mask", "vh_mask"): assert k in acq + def test_s3_uri_discovers_acquisitions(self) -> None: + """s3:// prefix is listed via s3fs; pathlib.glob is NOT used.""" + s3_files = [ + "bucket/prefix/s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC.tif", + "bucket/prefix/s1a_32TQM_vh_ASC_037_20230115t061234_GammaNaughtRTC.tif", + "bucket/prefix/s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC_BorderMask.tif", + "bucket/prefix/s1a_32TQM_vh_ASC_037_20230115t061234_GammaNaughtRTC_BorderMask.tif", + ] + with patch("s3fs.S3FileSystem.glob", return_value=s3_files): + acqs = discover_s1tiling_acquisitions("s3://bucket/prefix/") + assert len(acqs) == 1 + expected_vv = "s3://bucket/prefix/s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC.tif" + assert acqs[0]["vv"] == expected_vv + # ============================================================================= # Phase 3: Conditions ingestion tests @@ -788,3 +803,12 @@ def test_skips_non_matching(self, tmp_path: Path) -> None: def test_empty_directory(self, tmp_path: Path) -> None: conditions = discover_s1tiling_conditions(tmp_path) assert len(conditions) == 0 + + def test_s3_uri_discovers_conditions(self) -> None: + """s3:// prefix is listed via s3fs; pathlib.glob is NOT used.""" + s3_files = ["bucket/prefix/GAMMA_AREA_32TQM_037.tif"] + with patch("s3fs.S3FileSystem.glob", return_value=s3_files): + conditions = discover_s1tiling_conditions("s3://bucket/prefix/") + assert len(conditions) == 1 + assert conditions[0]["tile"] == "32TQM" + assert conditions[0]["orbit"] == "037" From 8d19118e73b9bf10f10b4094c65021dc46cc497f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Fri, 29 May 2026 11:39:06 +0100 Subject: [PATCH 03/19] fix(stac): point vv/vh asset hrefs to orbit group root; add tile_matrix_set Per Emmanuel's review on #173: - vv/vh asset hrefs now point to {store}/{orbit} (the zarr group that carries multiscales), matching how S2 reflectance assets point to measurements/reflectance. TiTiler reads tile_matrix_set from that group's multiscales attributes. - create_s1_store now writes tile_matrix_set into the orbit group's multiscales so TiTiler can discover the tiling scheme without error. - Updated test_asset_hrefs to assert the new orbit-group href. Co-Authored-By: Claude Sonnet 4.6 --- src/eopf_geozarr/conversion/s1_ingest.py | 44 ++++++++++++++++++++++++ src/eopf_geozarr/stac/s1_rtc.py | 4 +-- tests/test_s1_stac.py | 6 ++-- 3 files changed, 49 insertions(+), 5 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 938c6f40..6dab514c 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -22,6 +22,7 @@ import numpy as np import rasterio +import rasterio.crs import structlog import zarr import zarr.codecs @@ -296,6 +297,47 @@ def _create_spatial_coordinate_arrays( ) +def _create_tile_matrix_set(crs_string: str, bounds: list[float], layout: list[dict]) -> dict: + """Build an OGC TileMatrixSet for the orbit group from the store layout. + + Each resolution level becomes one tile matrix with matrixWidth=matrixHeight=1 + (one tile covers the full MGRS tile extent), matching the GeoZarr convention + used by S2 measurements/reflectance groups. + """ + left, bottom, right, top = bounds + native_crs = rasterio.crs.CRS.from_string(crs_string) + epsg = native_crs.to_epsg() + crs_uri = ( + f"http://www.opengis.net/def/crs/EPSG/0/{epsg}" + if epsg + else native_crs.to_wkt() + ) + tile_matrices = [] + for entry in layout: + h, w = entry["spatial:shape"] + cell_size = max((right - left) / w, (top - bottom) / h) + tile_matrices.append( + { + "id": entry["asset"], + "scaleDenominator": cell_size * 3779.5275, + "cellSize": cell_size, + "pointOfOrigin": [left, top], + "tileWidth": w, + "tileHeight": h, + "matrixWidth": 1, + "matrixHeight": 1, + } + ) + return { + "id": f"Native_CRS_{epsg or 'Custom'}", + "title": f"Native CRS Tile Matrix Set (EPSG:{epsg})", + "crs": crs_uri, + "supportedCRS": crs_uri, + "orderedAxes": ["X", "Y"], + "tileMatrices": tile_matrices, + } + + def create_s1_store( store_path: str | Path, orbit_direction: str, @@ -306,6 +348,7 @@ def create_s1_store( Returns the root group. """ layout = compute_multiscales_layout(metadata.shape, metadata.spatial_transform) + tile_matrix_set = _create_tile_matrix_set(metadata.crs, metadata.bounds, layout) root = zarr.open_group(str(store_path), mode="w-", zarr_format=3) orbit_group = root.create_group(orbit_direction) @@ -316,6 +359,7 @@ def create_s1_store( "multiscales": { "layout": layout, "resampling_method": "average", + "tile_matrix_set": tile_matrix_set, }, "proj:code": metadata.crs, "spatial:dimensions": ["y", "x"], diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index 3c57c567..79573f9c 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -137,10 +137,10 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: item.add_asset( pol, pystac.Asset( - href=f"{store_str}/{preferred_orbit}/r10m/{pol}", + href=f"{store_str}/{preferred_orbit}", media_type="application/vnd.zarr; version=3", roles=["data"], - title=f"{pol.upper()} polarisation (native 10 m)", + title=f"{pol.upper()} polarisation", ), ) diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index a67a9c87..513ce326 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -165,14 +165,14 @@ def test_empty_store_raises(tmp_path: Path) -> None: def test_asset_hrefs(tmp_path: Path) -> None: - """zarr-store href = store URI; vv/vh hrefs = {store}/{orbit}/r10m/{pol}.""" + """zarr-store href = store URI; vv/vh hrefs = {store}/{orbit} (orbit group root, per geozarr spec).""" store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") store_str = str(store) assert item.assets["zarr-store"].href == store_str - assert item.assets["vv"].href == f"{store_str}/descending/r10m/vv" - assert item.assets["vh"].href == f"{store_str}/descending/r10m/vh" + assert item.assets["vv"].href == f"{store_str}/descending" + assert item.assets["vh"].href == f"{store_str}/descending" def test_sar_extension_fields(tmp_path: Path) -> None: From e44c71cfd8063f54cd0496d1439fb3891efab722 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Fri, 29 May 2026 14:22:25 +0100 Subject: [PATCH 04/19] fix(ingest): add proj:code to resolution group attrs for TiTiler compatibility MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit TiTiler's _validate_zarr checks ds.rio.crs on the native resolution group (r10m). Without proj:code in the resolution group's attrs, rioxarray returns None and the group is excluded → groups=[] → bounds unpack error. S2 zarr already has proj:code at the resolution group level; mirror that for S1. Co-Authored-By: Claude Sonnet 4.6 --- src/eopf_geozarr/conversion/s1_ingest.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 6dab514c..442c84c1 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -376,6 +376,7 @@ def create_s1_store( { "spatial:shape": [level_h, level_w], "spatial:transform": level_entry["spatial:transform"], + "proj:code": metadata.crs, } ) From d1d1e31b8367562ab50219ced6a428877855ad72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Mon, 1 Jun 2026 10:38:23 +0100 Subject: [PATCH 05/19] fix(s1): write CF grid_mapping, drop tile_matrix_set from RTC store (#176) The deployed titiler-eopf reader resolves a group's CRS via rioxarray, which reads a CF spatial_ref/crs_wkt coordinate -- not the GeoZarr proj:code attr. S1 GRD RTC stores carried only proj:code, so rio.crs was None and every multiscale group was rejected (HTTP 500 on render/info). The S2 converter already writes a spatial_ref grid-mapping; this brings S1 in line. Also stop writing tile_matrix_set into the orbit-group multiscales: the data-model owner confirmed it is not part of the S1 GRD RTC data model, and the data_api reader already treats it as optional (MISSING default). - add _add_grid_mapping() using pyproj CRS.to_cf() (the same source rioxarray uses, so no hard-coded projection) and call it for every resolution level (both store-creation paths) and the conditions group - remove _create_tile_matrix_set() and its multiscales entry - tests: assert no tile_matrix_set, spatial_ref present + rio.crs resolves to the native EPSG, conditions group carries grid_mapping This fixes at the source what data-pipeline currently works around in ingest_v1_s1_rtc.py (_patch_cf_grid_mapping / _patch_tile_matrix_limits). Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/conversion/s1_ingest.py | 71 +++++++++++------------- tests/test_s1_rtc_ingest.py | 35 ++++++++++++ 2 files changed, 66 insertions(+), 40 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 442c84c1..7328ff4d 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -22,10 +22,10 @@ import numpy as np import rasterio -import rasterio.crs import structlog import zarr import zarr.codecs +from pyproj import CRS as PyprojCRS from zarr_cm import geo_proj from zarr_cm import multiscales as multiscales_cm from zarr_cm import spatial as spatial_cm @@ -297,45 +297,34 @@ def _create_spatial_coordinate_arrays( ) -def _create_tile_matrix_set(crs_string: str, bounds: list[float], layout: list[dict]) -> dict: - """Build an OGC TileMatrixSet for the orbit group from the store layout. +def _add_grid_mapping(group: zarr.Group, crs_string: str) -> None: + """Add a CF ``spatial_ref`` grid-mapping coordinate to a group holding (y, x) arrays. - Each resolution level becomes one tile matrix with matrixWidth=matrixHeight=1 - (one tile covers the full MGRS tile extent), matching the GeoZarr convention - used by S2 measurements/reflectance groups. + rioxarray -- and TiTiler's GeoZarr reader -- resolve the CRS from a CF + ``spatial_ref``/``crs_wkt`` coordinate; the GeoZarr ``proj:code`` attribute alone + is not read. This mirrors the S2 converter (which writes a ``spatial_ref`` + grid-mapping variable via ``rio.write_crs``). The CF attributes come from + ``pyproj.CRS.to_cf()`` -- the same source rioxarray uses -- so the projection is + described correctly for any CRS rather than hard-coded. + + The scalar ``spatial_ref`` array is created (if absent) and every (y, x) data + array in the group is given ``grid_mapping = "spatial_ref"``. """ - left, bottom, right, top = bounds - native_crs = rasterio.crs.CRS.from_string(crs_string) - epsg = native_crs.to_epsg() - crs_uri = ( - f"http://www.opengis.net/def/crs/EPSG/0/{epsg}" - if epsg - else native_crs.to_wkt() - ) - tile_matrices = [] - for entry in layout: - h, w = entry["spatial:shape"] - cell_size = max((right - left) / w, (top - bottom) / h) - tile_matrices.append( - { - "id": entry["asset"], - "scaleDenominator": cell_size * 3779.5275, - "cellSize": cell_size, - "pointOfOrigin": [left, top], - "tileWidth": w, - "tileHeight": h, - "matrixWidth": 1, - "matrixHeight": 1, - } - ) - return { - "id": f"Native_CRS_{epsg or 'Custom'}", - "title": f"Native CRS Tile Matrix Set (EPSG:{epsg})", - "crs": crs_uri, - "supportedCRS": crs_uri, - "orderedAxes": ["X", "Y"], - "tileMatrices": tile_matrices, - } + cf_attrs = PyprojCRS.from_user_input(crs_string).to_cf() + # rioxarray writes both ``crs_wkt`` and a ``spatial_ref`` attr holding the WKT. + cf_attrs["spatial_ref"] = cf_attrs["crs_wkt"] + cf_attrs["_ARRAY_DIMENSIONS"] = [] + + if "spatial_ref" in group: + sref = group["spatial_ref"] + else: + sref = group.create_array("spatial_ref", shape=(), dtype="int64", fill_value=0) + sref[...] = 0 + sref.attrs.update(cf_attrs) + + for name, arr in group.arrays(): + if name != "spatial_ref" and {"y", "x"}.issubset(arr.metadata.dimension_names or ()): + arr.attrs["grid_mapping"] = "spatial_ref" def create_s1_store( @@ -348,7 +337,6 @@ def create_s1_store( Returns the root group. """ layout = compute_multiscales_layout(metadata.shape, metadata.spatial_transform) - tile_matrix_set = _create_tile_matrix_set(metadata.crs, metadata.bounds, layout) root = zarr.open_group(str(store_path), mode="w-", zarr_format=3) orbit_group = root.create_group(orbit_direction) @@ -359,7 +347,6 @@ def create_s1_store( "multiscales": { "layout": layout, "resampling_method": "average", - "tile_matrix_set": tile_matrix_set, }, "proj:code": metadata.crs, "spatial:dimensions": ["y", "x"], @@ -406,6 +393,7 @@ def create_s1_store( _create_spatial_coordinate_arrays( level_group, level_h, level_w, level_entry["spatial:transform"] ) + _add_grid_mapping(level_group, metadata.crs) # Coordinate variables at native resolution only r10m = orbit_group["r10m"] @@ -588,6 +576,7 @@ def ingest_s1tiling_acquisition( _create_spatial_coordinate_arrays( level_group, level_h, level_w, level_entry["spatial:transform"] ) + _add_grid_mapping(level_group, meta.crs) r10m = orbit_group["r10m"] for name, dtype, fill in [ ("time", "int64", 0), @@ -981,6 +970,8 @@ def ingest_s1tiling_conditions( max=float(np.nanmax(data)), ) + _add_grid_mapping(conditions, ref_crs) + log.info( "Conditions ingestion complete", orbit_direction=orbit_direction, diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py index e404d642..0414dbff 100644 --- a/tests/test_s1_rtc_ingest.py +++ b/tests/test_s1_rtc_ingest.py @@ -244,6 +244,38 @@ def test_array_metadata(self, s1_store_path: Path, sample_metadata: S1TilingMeta assert r10m["vv"].dtype == np.float32 assert r10m["border_mask"].dtype == np.uint8 + def test_no_tile_matrix_set( + self, s1_store_path: Path, sample_metadata: S1TilingMetadata + ) -> None: + # tile_matrix_set is not part of the S1 GRD RTC data model (confirmed with the + # data-model owner): the multiscales attribute must not carry one. + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + ms = dict(root["ascending"].attrs)["multiscales"] + assert "tile_matrix_set" not in ms + assert "layout" in ms + + def test_cf_grid_mapping_resolves_crs( + self, s1_store_path: Path, sample_metadata: S1TilingMetadata + ) -> None: + # Each resolution level carries a CF spatial_ref grid-mapping so rioxarray (and + # TiTiler's GeoZarr reader) can resolve the CRS -- the geozarr proj:code attr + # alone is not read by rioxarray. + import rioxarray # noqa: F401 -- registers the .rio accessor + + create_s1_store(s1_store_path, "ascending", sample_metadata) + r10m = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3)["ascending"]["r10m"] + assert "spatial_ref" in list(r10m.array_keys()) + assert dict(r10m["vv"].attrs).get("grid_mapping") == "spatial_ref" + assert dict(r10m["vh"].attrs).get("grid_mapping") == "spatial_ref" + + ds = xr.open_zarr( + str(s1_store_path / "ascending" / "r10m"), + consolidated=False, + decode_coords="all", + ) + assert ds.rio.crs is not None + assert ds.rio.crs.to_epsg() == 32633 + def test_coordinate_variables( self, s1_store_path: Path, sample_metadata: S1TilingMetadata ) -> None: @@ -609,6 +641,9 @@ def test_conditions_group_attributes( assert attrs["spatial:dimensions"] == ["y", "x"] assert len(attrs["spatial:transform"]) == 6 assert attrs["spatial:shape"] == [SIZE, SIZE] + # CF grid-mapping so rioxarray can resolve the CRS of the condition arrays + assert "spatial_ref" in list(conditions.array_keys()) + assert dict(conditions["gamma_area_037"].attrs).get("grid_mapping") == "spatial_ref" def test_gamma_area_array_shape_and_dtype( self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path From 86206dd7ef3aa310dd803bb186789afe690af8b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Tue, 9 Jun 2026 13:57:36 +0100 Subject: [PATCH 06/19] feat(s1-stac): add render extension with dual-pol RGB composite preview MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous S1 RTC preview rendered a single VH band as grayscale with an incorrect rescale (0,219), unsuitable for linear gamma0 RTC values. Declare a render-extension `renders.rgb` config producing the standard dual-pol false- colour composite (R=VV, G=VH, B=VV/VH ratio) with rescale 0–0.1 and tilesize 256, referencing the preferred orbit group. Downstream titiler-based services consume this to generate previews/tiles/tilejson. Co-Authored-By: Claude Opus 4.8 --- src/eopf_geozarr/stac/s1_rtc.py | 23 ++++++++++++++++++++++- tests/test_s1_stac.py | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 1 deletion(-) diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index 79573f9c..23b8a0f8 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -14,10 +14,29 @@ SAR_EXT = "https://stac-extensions.github.io/sar/v1.0.0/schema.json" SAT_EXT = "https://stac-extensions.github.io/sat/v1.0.0/schema.json" PROJ_EXT = "https://stac-extensions.github.io/projection/v2.0.0/schema.json" +RENDER_EXT = "https://stac-extensions.github.io/render/v1.0.0/schema.json" _ORBIT_PREFERENCE = ("ascending", "descending") +def _rgb_render(orbit: str) -> dict[str, object]: + """Build the dual-pol RGB composite render config for the given orbit group. + + Produces a 3-band false-colour composite (R=VV, G=VH, B=VV/VH ratio) that + titiler renders into previews/tiles. ``bidx=[1]`` selects the single time + slice from each multi-band variable; ``rescale`` is in linear gamma0 units. + """ + vv = f"/{orbit}:vv" + vh = f"/{orbit}:vh" + return { + "title": "VV, VH, VV/VH composite", + "expression": f"{vv};{vh};({vv})/({vh})", + "rescale": [[0.0, 0.1], [0.0, 0.1], [0.0, 0.1]], + "bidx": [1], + "tilesize": 256, + } + + def _utm_to_wgs84( proj_code: str, utm_bbox: list[float] ) -> tuple[float, float, float, float]: @@ -118,8 +137,10 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: "sat:orbit_state": preferred_orbit, # Projection extension "proj:code": preferred_proj_code, + # Render extension: dual-pol RGB composite for previews/tiles + "renders": {"rgb": _rgb_render(preferred_orbit)}, }, - stac_extensions=[SAR_EXT, SAT_EXT, PROJ_EXT], + stac_extensions=[SAR_EXT, SAT_EXT, PROJ_EXT, RENDER_EXT], collection=collection_id, ) diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index 513ce326..2a40eb06 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -189,3 +189,36 @@ def test_sar_extension_fields(tmp_path: Path) -> None: sar_ext_uri = "https://stac-extensions.github.io/sar/v1.0.0/schema.json" assert sar_ext_uri in item.stac_extensions + + +def test_render_extension_rgb_composite(tmp_path: Path) -> None: + """Item must declare a render-extension RGB composite using the preferred orbit.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + render_ext_uri = "https://stac-extensions.github.io/render/v1.0.0/schema.json" + assert render_ext_uri in item.stac_extensions + + rgb = item.properties["renders"]["rgb"] + assert rgb["expression"] == "/descending:vv;/descending:vh;(/descending:vv)/(/descending:vh)" + assert rgb["rescale"] == [[0.0, 0.1], [0.0, 0.1], [0.0, 0.1]] + assert rgb["bidx"] == [1] + assert rgb["tilesize"] == 256 + + +def test_render_uses_ascending_when_preferred(tmp_path: Path) -> None: + """When ascending is the preferred orbit, the render expression must reference it.""" + store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + for orbit_dir in ("descending", "ascending"): + og = root.create_group(orbit_dir) + og.attrs.update({"proj:code": CRS, "spatial:bbox": UTM_BBOX}) + r10m = og.create_group("r10m") + t_arr = r10m.create_array("time", shape=(1,), dtype="int64", chunks=(512,)) + t_arr[:] = [T1_NS] + p_arr = r10m.create_array("platform", shape=(1,), dtype=" Date: Tue, 9 Jun 2026 14:31:16 +0100 Subject: [PATCH 07/19] refactor(s1-stac): emit a single rescale pair in the RGB render A single [min, max] pair applies to all bands in titiler, so emitting three identical pairs was redundant (and forced collapse logic in the consumer). Same rendered output. Co-Authored-By: Claude Opus 4.8 --- src/eopf_geozarr/stac/s1_rtc.py | 5 +++-- tests/test_s1_stac.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index 23b8a0f8..027ec423 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -24,14 +24,15 @@ def _rgb_render(orbit: str) -> dict[str, object]: Produces a 3-band false-colour composite (R=VV, G=VH, B=VV/VH ratio) that titiler renders into previews/tiles. ``bidx=[1]`` selects the single time - slice from each multi-band variable; ``rescale`` is in linear gamma0 units. + slice from each multi-band variable; the single ``rescale`` pair (linear + gamma0 units) is applied to every band. """ vv = f"/{orbit}:vv" vh = f"/{orbit}:vh" return { "title": "VV, VH, VV/VH composite", "expression": f"{vv};{vh};({vv})/({vh})", - "rescale": [[0.0, 0.1], [0.0, 0.1], [0.0, 0.1]], + "rescale": [[0.0, 0.1]], "bidx": [1], "tilesize": 256, } diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index 2a40eb06..9fcaa990 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -201,7 +201,7 @@ def test_render_extension_rgb_composite(tmp_path: Path) -> None: rgb = item.properties["renders"]["rgb"] assert rgb["expression"] == "/descending:vv;/descending:vh;(/descending:vv)/(/descending:vh)" - assert rgb["rescale"] == [[0.0, 0.1], [0.0, 0.1], [0.0, 0.1]] + assert rgb["rescale"] == [[0.0, 0.1]] assert rgb["bidx"] == [1] assert rgb["tilesize"] == 256 From 4c47cd5f9fa4544ab586e7db544b000ced22dc0d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Mon, 15 Jun 2026 10:00:29 +0100 Subject: [PATCH 08/19] fix(s1-rtc): TEMPORARY name store as s1-rtc-{tile} to match titiler render path (#246) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit titiler-eopf reconstructs the store path as s3://{bucket}/tests-output/{collection}/{item_id}.zarr and ignores the STAC asset href, so the GeoZarr cube must be named after the item id (s1-rtc-{tile}) for new tiles to preview. Parse the tile from the s1-rtc- prefix (was s1-grd-rtc-); the item id (s1-rtc-{tile}) is unchanged, so filename == item id by construction. Pairs with the data-pipeline direct-write change (store written at the tests-output path) which replaces the #250 auto-copy. Dev-phase, no consumers — disposable test data. Revert to s1-grd-rtc- when titiler-eopf#108 (resolve store from href) lands. Co-Authored-By: Claude Opus 4.8 --- src/eopf_geozarr/stac/s1_rtc.py | 5 ++++- tests/test_s1_stac.py | 12 +++++++----- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index 027ec423..297a6c20 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -69,7 +69,10 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: ValueError If the store contains no acquisitions. """ - tile_id = Path(zarr_store).name.removeprefix("s1-grd-rtc-").removesuffix(".zarr") + # TEMPORARY (#246): the store is written as s1-rtc-{tile}.zarr so its filename equals + # the item id, which titiler-eopf reconstructs as the render path (it ignores the asset + # href). Revert this prefix to "s1-grd-rtc-" when titiler-eopf#108 lands. + tile_id = Path(zarr_store).name.removeprefix("s1-rtc-").removesuffix(".zarr") root = zarr.open_consolidated(zarr_store, zarr_format=3) diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index 9fcaa990..fdfe5336 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -45,7 +45,9 @@ def _make_s1_store( """ if utm_bbox is None: utm_bbox = UTM_BBOX - store_path = tmp_path / f"s1-grd-rtc-{tile_id}.zarr" + # TEMPORARY (#246): store basename == item id (s1-rtc-{tile}) so titiler's reconstructed + # render path resolves; revert to "s1-grd-rtc-" when titiler-eopf#108 lands. + store_path = tmp_path / f"s1-rtc-{tile_id}.zarr" root = zarr.open_group(str(store_path), mode="w", zarr_format=3) for orbit_dir, acquisitions in orbits.items(): og = root.create_group(orbit_dir) @@ -105,7 +107,7 @@ def test_bbox_wgs84_from_utm(tmp_path: Path) -> None: def test_both_orbits_bbox_union(tmp_path: Path) -> None: """When ascending and descending are both present, the WGS84 bbox is the union.""" # Give ascending a different UTM bbox (shifted east) - store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + store_path = tmp_path / "s1-rtc-31TCH.zarr" root = zarr.open_group(str(store_path), mode="w", zarr_format=3) for orbit_dir, bbox in [ @@ -131,7 +133,7 @@ def test_both_orbits_bbox_union(tmp_path: Path) -> None: def test_ascending_preferred_for_assets(tmp_path: Path) -> None: """When both orbits present, vv/vh assets must point to the ascending group.""" - store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + store_path = tmp_path / "s1-rtc-31TCH.zarr" root = zarr.open_group(str(store_path), mode="w", zarr_format=3) for orbit_dir in ("descending", "ascending"): og = root.create_group(orbit_dir) @@ -150,7 +152,7 @@ def test_ascending_preferred_for_assets(tmp_path: Path) -> None: def test_empty_store_raises(tmp_path: Path) -> None: """A store with an orbit group but no acquisitions must raise ValueError.""" - store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + store_path = tmp_path / "s1-rtc-31TCH.zarr" root = zarr.open_group(str(store_path), mode="w", zarr_format=3) og = root.create_group("descending") og.attrs.update({"proj:code": CRS, "spatial:bbox": UTM_BBOX}) @@ -208,7 +210,7 @@ def test_render_extension_rgb_composite(tmp_path: Path) -> None: def test_render_uses_ascending_when_preferred(tmp_path: Path) -> None: """When ascending is the preferred orbit, the render expression must reference it.""" - store_path = tmp_path / "s1-grd-rtc-31TCH.zarr" + store_path = tmp_path / "s1-rtc-31TCH.zarr" root = zarr.open_group(str(store_path), mode="w", zarr_format=3) for orbit_dir in ("descending", "ascending"): og = root.create_group(orbit_dir) From b839a39c8dc5412d6511d94d39b73897a6e83fc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Mon, 15 Jun 2026 21:43:49 +0100 Subject: [PATCH 09/19] fix(s1-rtc): build STAC item from a non-consolidated cube store MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit build_s1_rtc_stac_item opened the cube with zarr.open_consolidated, which raises `ValueError: Consolidated metadata ... not found` when the store lacks root consolidated metadata. A per-tile cube grown by appending a time-slice to an existing same-orbit group ends up in exactly that state — re-consolidating an append on the S3 store is unreliable — so live STAC registration failed on the 2nd+ same-orbit acquisition for a tile (it broke the Pyrenees S1 RTC soak). The builder must not require consolidated metadata: titiler reads these stores fine without it. Fall back to a direct zarr.open_group(mode="r") read when the consolidated metadata is absent; re-raise any other ValueError. Adds a regression test: build from a 2-slice same-orbit store created with consolidate=False (fails before the fix in open_consolidated, passes after). Co-Authored-By: Claude Opus 4.8 --- src/eopf_geozarr/stac/s1_rtc.py | 11 ++++++++++- tests/test_s1_stac.py | 25 +++++++++++++++++++++++-- 2 files changed, 33 insertions(+), 3 deletions(-) diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index 297a6c20..eba7da6a 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -74,7 +74,16 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: # href). Revert this prefix to "s1-grd-rtc-" when titiler-eopf#108 lands. tile_id = Path(zarr_store).name.removeprefix("s1-rtc-").removesuffix(".zarr") - root = zarr.open_consolidated(zarr_store, zarr_format=3) + # A cube grown by appending a time-slice to an *existing same-orbit* group can end up without + # root consolidated metadata (re-consolidating an append on the S3 store is unreliable). The + # builder must not require it — fall back to reading the hierarchy directly, exactly as titiler + # does. See the data-model issue on the S1 RTC consolidated-metadata regression. + try: + root = zarr.open_consolidated(zarr_store, zarr_format=3) + except ValueError as exc: + if "consolidated metadata" not in str(exc).lower(): + raise + root = zarr.open_group(zarr_store, mode="r", zarr_format=3) all_times_ns: list[int] = [] wgs84_bboxes: list[tuple[float, float, float, float]] = [] diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index fdfe5336..cb2b9226 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -37,11 +37,14 @@ def _make_s1_store( tile_id: str = "31TCH", crs: str = CRS, utm_bbox: list[float] | None = None, + consolidate: bool = True, ) -> Path: - """Create a minimal consolidated S1 Zarr store. + """Create a minimal S1 Zarr store. ``orbits`` maps orbit_direction -> list of (time_ns, platform) tuples. Creates only the attrs and coordinate arrays needed by build_s1_rtc_stac_item. + ``consolidate=False`` skips writing root consolidated metadata, mirroring a cube that grew by + appending a time-slice to an existing same-orbit group (the builder must still read it). """ if utm_bbox is None: utm_bbox = UTM_BBOX @@ -59,7 +62,8 @@ def _make_s1_store( t_arr[:] = times p_arr = r10m.create_array("platform", shape=platforms.shape, dtype=" None: assert item.id == "s1-rtc-31TCH" +def test_builds_from_non_consolidated_store(tmp_path: Path) -> None: + """Regression: the builder must read a store that lacks root consolidated metadata. + + A per-tile cube grown by appending a time-slice to an existing same-orbit group can end up + without root consolidated metadata (re-consolidating an S3 append is unreliable), which made + ``zarr.open_consolidated`` raise ``ValueError: Consolidated metadata ... not found`` and broke + STAC registration in the live S1 RTC pipeline. The builder must fall back to a direct read. + """ + store = _make_s1_store( + tmp_path, {"ascending": [(T1_NS, "S1A"), (T2_NS, "S1A")]}, consolidate=False + ) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + assert item.id == "s1-rtc-31TCH" + assert item.properties["start_datetime"] + assert item.properties["end_datetime"] + + def test_temporal_range_min_max(tmp_path: Path) -> None: """start_datetime/end_datetime must span the full time range across all acquisitions.""" store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A"), (T2_NS, "S1A")]}) From 12a3df71183c1122ad51610c77f02b4768bd381b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Wed, 17 Jun 2026 12:29:30 +0100 Subject: [PATCH 10/19] fix(s1-rtc): CF-encode `time` at every level for datetime-based rendering (#192) (#193) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Per-acquisition previews render via titiler `sel=time={INTEGER index}` — a positional index that breaks once a cube's time axis goes non-monotonic (a cross-run append of an earlier-dated scene; proved on 31TEH: `[06-08, 06-07]` → the 06-07 item has no preview). Keeping positional indices correct would need cube reorder + full re-registration on every append (not scalable). The deployed titiler (v0.5.0) already does label-based `da.sel(..., method=...)` via `open_datatree(decode_times=True)`; the only gap is the cube's `time` array, which is a bare int64 with no CF datetime metadata, so it can't be a datetime index. Encoding it lets previews select by exact datetime, which works even on a non-monotonic axis — no reorder, register only the new item. titiler picks a multiscale level by zoom (previews use a coarse level), so `time` must resolve there. Putting it only at the orbit-group level and inheriting fails to open while r10m stays bare int64 (AlignmentError: int64 vs datetime64 on the shared dim). So encode `time` consistently at EVERY level: - new `_create_time_coordinate_array` + `TIME_CF_ATTRS` (units/calendar/ standard_name); created at every level in create + append-group paths. - the append write resizes/writes `time` at every level (was r10m only). - absolute_orbit/relative_orbit/platform stay at r10m (not selected on). - stored dtype stays int64 ns, so register_per_acquisition's raw read is unaffected. Adds TestTimeCFDatetime incl. exact datetime `.sel` on a non-monotonic axis at the native and a coarse level. 55 passed. Refs EOPF-Explorer/data-model#192 Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/conversion/s1_ingest.py | 54 +++++++++++--- tests/test_s1_rtc_ingest.py | 95 ++++++++++++++++++++++++ 2 files changed, 140 insertions(+), 9 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 7328ff4d..751b0913 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -297,6 +297,37 @@ def _create_spatial_coordinate_arrays( ) +# CF datetime encoding for the `time` coordinate. Without it the array is a bare int64 and readers +# (xarray / TiTiler's `open_datatree(decode_times=True)`) cannot expose `time` as a datetime index, so +# per-acquisition previews can only select positionally (`sel=time={index}`) — fragile once a cube's +# time axis goes non-monotonic. With these attrs `time` decodes to datetime64 and previews can render by +# `sel=time={datetime}` (order-immune). The stored dtype stays int64 nanoseconds. See data-model #192. +TIME_CF_ATTRS = { + "units": "nanoseconds since 1970-01-01", + "calendar": "proleptic_gregorian", + "standard_name": "time", +} + + +def _create_time_coordinate_array(level_group: zarr.Group) -> None: + """Create the CF-encoded ``time`` coordinate (length 0, grown on append) on one level group. + + Replicated at EVERY multiscale level (with identical values) so datetime ``.sel`` resolves at + whatever level a reader renders — TiTiler picks a coarse level for previews, and a level lacking a + ``time`` coordinate cannot be selected by datetime. Keeping the dtype/values consistent across + levels is also required for the datatree to open (mixed int64/datetime64 fails alignment). + """ + t_arr = level_group.create_array( + "time", + shape=(0,), + dtype="int64", + chunks=(512,), + fill_value=0, + dimension_names=["time"], + ) + t_arr.attrs.update({**TIME_CF_ATTRS, "_ARRAY_DIMENSIONS": ["time"]}) + + def _add_grid_mapping(group: zarr.Group, crs_string: str) -> None: """Add a CF ``spatial_ref`` grid-mapping coordinate to a group holding (y, x) arrays. @@ -394,11 +425,12 @@ def create_s1_store( level_group, level_h, level_w, level_entry["spatial:transform"] ) _add_grid_mapping(level_group, metadata.crs) + # `time` coordinate on every level so datetime `.sel` resolves at any rendered scale (#192). + _create_time_coordinate_array(level_group) - # Coordinate variables at native resolution only + # Per-time metadata coordinates at native resolution only (not selected on by readers). r10m = orbit_group["r10m"] for name, dtype, fill in [ - ("time", "int64", 0), ("absolute_orbit", "int32", 0), ("relative_orbit", "int32", 0), ]: @@ -577,9 +609,10 @@ def ingest_s1tiling_acquisition( level_group, level_h, level_w, level_entry["spatial:transform"] ) _add_grid_mapping(level_group, meta.crs) + # `time` coordinate on every level (datetime `.sel` at any scale, #192). + _create_time_coordinate_array(level_group) r10m = orbit_group["r10m"] for name, dtype, fill in [ - ("time", "int64", 0), ("absolute_orbit", "int32", 0), ("relative_orbit", "int32", 0), ]: @@ -649,7 +682,10 @@ def ingest_s1tiling_acquisition( log.info("Overviews generated", levels=len(data_by_level)) - # Write data at all levels + dt_ns = np.datetime64(meta.datetime).astype("datetime64[ns]").astype(np.int64) + + # Write data + the `time` coordinate at all levels (time is replicated per level so datetime + # `.sel` resolves at any rendered scale, #192). for level_name, (vv_lev, vh_lev, mask_lev) in data_by_level.items(): level = orbit[level_name] h, w = vv_lev.shape @@ -662,12 +698,12 @@ def ingest_s1tiling_acquisition( level["vh"][current_size, :, :] = vh_lev level["border_mask"][current_size, :, :] = mask_lev - # Append coordinate variables - for coord_name in ["time", "absolute_orbit", "relative_orbit", "platform"]: - r10m[coord_name].resize((new_size,)) + level["time"].resize((new_size,)) + level["time"][current_size] = dt_ns - dt_ns = np.datetime64(meta.datetime).astype("datetime64[ns]").astype(np.int64) - r10m["time"][current_size] = dt_ns + # Per-time metadata coordinates at native resolution only. + for coord_name in ["absolute_orbit", "relative_orbit", "platform"]: + r10m[coord_name].resize((new_size,)) r10m["absolute_orbit"][current_size] = meta.absolute_orbit r10m["relative_orbit"][current_size] = meta.relative_orbit r10m["platform"][current_size] = meta.platform diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py index 0414dbff..7fea3229 100644 --- a/tests/test_s1_rtc_ingest.py +++ b/tests/test_s1_rtc_ingest.py @@ -847,3 +847,98 @@ def test_s3_uri_discovers_conditions(self) -> None: assert len(conditions) == 1 assert conditions[0]["tile"] == "32TQM" assert conditions[0]["orbit"] == "037" + + +# ============================================================================= +# CF datetime `time` coordinate — render-by-datetime support (data-model #192) +# ============================================================================= + +_LEVELS = ["r10m", "r20m", "r60m", "r120m", "r360m", "r720m"] + + +class TestTimeCFDatetime: + """`time` is CF-encoded at every multiscale level so readers decode it to datetime64 and can + select a slice by datetime (`sel=time={datetime}`) at any rendered scale, even on a non-monotonic + axis. This replaces the fragile positional `sel=time={index}` rendering (#192).""" + + def _paths(self, geotiff_dir: Path, stamp: str) -> tuple[Path, Path, Path]: + vv = geotiff_dir / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC.tif" + vh = geotiff_dir / f"s1a_32TQM_vh_ASC_037_{stamp}_GammaNaughtRTC.tif" + mask = geotiff_dir / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC_BorderMask.tif" + return vv, vh, mask + + def test_time_has_cf_attrs_at_every_level( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + vv, vh, mask = self._paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + for level in _LEVELS: + attrs = dict(root["ascending"][level]["time"].attrs) + assert attrs.get("units") == "nanoseconds since 1970-01-01", level + assert attrs.get("calendar") == "proleptic_gregorian", level + assert attrs.get("standard_name") == "time", level + + def test_open_datatree_decodes_time_to_datetime64( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + vv, vh, mask = self._paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") + dt = xr.open_datatree( + str(s1_store_path), engine="zarr", decode_times=True, consolidated=False + ) + for level in ("r10m", "r720m"): + da = dt["ascending"][level]["vv"] + assert "time" in da.coords, level + assert np.issubdtype(da["time"].dtype, np.datetime64), level + + def test_exact_datetime_sel_on_nonmonotonic_axis( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """Ingest LATER acq first then EARLIER → non-monotonic axis; exact datetime `.sel` still + returns the right physical slice at both the native and a coarse level (the 31TEH case).""" + vv2, vh2, mask2 = self._paths(s1_geotiff_dir, "20230127t061235") # 2023-01-27 (later) + vv1, vh1, mask1 = self._paths(s1_geotiff_dir, "20230115t061234") # 2023-01-15 (earlier) + ingest_s1tiling_acquisition(vv2, vh2, mask2, s1_store_path, "ascending") # -> index 0 + ingest_s1tiling_acquisition(vv1, vh1, mask1, s1_store_path, "ascending") # -> index 1 + + dt = xr.open_datatree( + str(s1_store_path), engine="zarr", decode_times=True, consolidated=False + ) + times = dt["ascending"]["r10m"]["time"].values + assert times[0] > times[1], "axis should be non-monotonic (later acq appended first)" + + early = np.datetime64("2023-01-15T06:12:34") # physical index 1 + later = np.datetime64("2023-01-27T06:12:35") # physical index 0 + for level in ("r10m", "r720m"): + vvda = dt["ascending"][level]["vv"] + np.testing.assert_array_equal( + vvda.sel(time=early).values, vvda.isel(time=1).values, err_msg=f"{level} early" + ) + np.testing.assert_array_equal( + vvda.sel(time=later).values, vvda.isel(time=0).values, err_msg=f"{level} later" + ) + + def test_time_values_identical_across_levels( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + vv1, vh1, mask1 = self._paths(s1_geotiff_dir, "20230115t061234") + vv2, vh2, mask2 = self._paths(s1_geotiff_dir, "20230127t061235") + ingest_s1tiling_acquisition(vv1, vh1, mask1, s1_store_path, "ascending") + ingest_s1tiling_acquisition(vv2, vh2, mask2, s1_store_path, "ascending") + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + ref = np.asarray(root["ascending"]["r10m"]["time"][:]) + assert ref.shape == (2,) + for level in _LEVELS[1:]: + np.testing.assert_array_equal(np.asarray(root["ascending"][level]["time"][:]), ref) + + def test_r10m_time_still_int64_for_register( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """register_per_acquisition reads r10m/time as raw int64 ns — CF attrs must not change that.""" + vv, vh, mask = self._paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + arr = root["ascending"]["r10m"]["time"] + assert arr.dtype == np.dtype("int64") + assert str(np.datetime64(int(arr[0]), "ns")).startswith("2023-01-15") From 930131c22f5bd91600a956ca4fffe772504cee7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Thu, 18 Jun 2026 16:45:23 +0100 Subject: [PATCH 11/19] fix(s1-rtc): enrich STAC items & consolidate construction into the library (#195) (#196) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * fix(s1-rtc): enrich cube STAC item + fix vv/vh asset model (#195) The S1 RTC cube STAC item carried duplicate vv/vh assets with byte-identical hrefs (bug #1), exposed no asset for the descending orbit group of a dual-orbit cube (bug #2), and was missing most of the metadata the S2 L2A reference carries. Rework `build_s1_rtc_stac_item`: - Asset model: replace the two indistinguishable vv/vh assets with one `gamma0-rtc-backscatter-{asc,desc}` asset per present orbit group, exposing VV/VH as named STAC 1.1 `bands` (no more duplicate hrefs), plus a `border-mask-{orbit}` asset for the valid-data mask. Each γ⁰ asset carries data_type/nodata/unit/gsd. - Identity: add `constellation`, `instruments: [c-sar]`, `gsd` (platform stays per-acquisition — a cube can mix S1A/S1C). - Projection: add proj:bbox, and proj:shape/proj:transform (best-effort from the r10m group attrs). - Datacube extension: cube:dimensions (time as a bounded extent — not a values list, the cube grows by appending — plus x/y) + cube:variables. - Descriptive: title/description; created (earliest acquisition, stable across rebuilds) / updated (build time) via the timestamps extension. - Fold the corrected render rescale (0.0,0.2) into `_rgb_render`. Out of scope (tracked follow-ups): statistics (#157), processing:software / DEM lineage, and per-product CDSE fields are not in the store; per-acquisition item construction is moved in a follow-up commit. Co-Authored-By: Claude Opus 4.8 * feat(s1-rtc): move per-acquisition + coverage STAC construction into the library Consolidate the *pure* S1 RTC STAC construction in eopf_geozarr.stac.s1_rtc so item shape has a single owner (the library); registration (STAC-API upsert, S3 alternates, gateway rewrite, TiTiler links, CDSE discovery) stays in data-pipeline. Adds, ported from the data-pipeline scripts (deployment URLs stripped): - acquisition_id(tile_id, when) - build_s1_rtc_per_acquisition_items(store, *, orbit, collection_id) -> list[pystac.Item]: one single-datetime item per cube time slice of an orbit, reoriented to that orbit (sat:orbit_state + render expression + only that orbit's γ⁰/mask assets), per-slice platform, datacube structure dropped (a single acquisition is not a cube). Deployment-agnostic: carries the render config + datetime; the registration layer derives the cube-endpoint TiTiler links (sel=time={datetime}) from it. - Slice / pick_slice / slice_coverages: coverage-based preview-slice selection (reads border_mask at r720m). Ported the construction unit tests (pick_slice, slice_coverages, and the per-acquisition construction subset); the link/thumbnail/alternate-asset and rescale-override tests stay in data-pipeline as registration concerns. Follow-up (separate data-pipeline PR): bump the data-model pin, delete the moved functions from register_per_acquisition.py / register_v1_s1_rtc.py, and rewire them to consume this library output. Co-Authored-By: Claude Opus 4.8 * fix(s1-rtc): per-acquisition items carry their own orbit footprint Review follow-up. A per-acquisition item cloned the cube base, so on a dual-orbit cube every item inherited the WGS84 union bbox/geometry of both orbits and the preferred (ascending) orbit's proj:bbox — advertising a footprint wider than, and a data bbox from the wrong orbit relative to, the acquisition's own orbit. Recompute bbox/geometry/proj:bbox from the run orbit's spatial:bbox (factored the bbox->polygon construction into `_bbox_to_geometry`, shared with the cube builder), and give per-acquisition items a single-acquisition description rather than the cube's "datacube" wording. Tests: per-acq footprint = run orbit (not the union); `created` = earliest acquisition (idempotent across rebuilds); ValueError on an invalid or absent orbit. Co-Authored-By: Claude Opus 4.8 * fix(s1-rtc): omit orbit_state on dual-orbit cubes; expose datacube dimension sizes Two review observations from inspecting the live items: - sat:orbit_state is single-valued, so on a dual-orbit cube (both ascending and descending groups present) a single value mislabels half the slices. Set it (and declare the SAT extension) only for a single-orbit cube; per-acquisition items, which are single-orbit, still carry the real per-orbit value. Ensure per-acq items declare the SAT extension even when cloned from a dual-orbit base that omitted it. - cube:dimensions previously carried only `extent`, so the number of elements per dimension wasn't visible. The time axis is irregularly sampled, so list its discrete `values` (count = number of acquisitions); the regular x/y axes get a `step` (size derivable from extent; exact pixel count is also in proj:shape). Co-Authored-By: Claude Opus 4.8 * fix(s1-rtc): datacube time axis carries extent only (drop values enumeration) Enumerating every acquisition as cube:dimensions.time.values does not scale as the cube is appended over the mission. No STAC extension provides a scalable temporal element-count field (datacube temporal_dimension offers only values/step/extent, and step is null for S1's irregular sampling; projection proj:shape is spatial-only; raster has no dimension shape). So the time axis carries a bounded `extent` only. x/y element counts remain available via proj:shape (exact) and the datacube x/y extent+step. The number of acquisitions is obtained from the per-acquisition collection (STAC search count), which stays correct as the cube grows. Co-Authored-By: Claude Opus 4.8 * fix(s1-rtc): correct datacube variable_type, platform value, created, time values Review of the live items surfaced several metadata issues: - cube:variables used `type` (a non-standard, silently-ignored key) instead of the datacube extension's `variable_type`; border_mask is now `auxiliary`, vv/vh `data`. - per-acquisition `platform` carried the store's short code (e.g. `s1a`); normalize to the STAC convention `sentinel-1a` (mirrors the Sentinel-2 reference). - `created` (timestamps ext) was set to an acquisition time, which misuses the field (it means the item's creation instant); omit it and keep `updated` = build time. - the cube's irregular time axis now lists its discrete `values` (the acquisition instants) so the number of time steps is visible; the regular x/y axes stay extent + step (their pixel count is in proj:shape). Co-Authored-By: Claude Opus 4.8 * refactor(s1-rtc): robust EPSG parse for datacube reference_system + clarify per-acq proj inheritance Review follow-ups (no behaviour change for the EPSG:NNNN stores in use): - derive the datacube reference_system via pyproj.CRS(...).to_epsg() instead of splitting the proj:code string, so non-"EPSG:NNNN" CRS forms don't break it. - comment why per-acquisition items recompute proj:bbox from their orbit but inherit proj:code/shape/transform (the MGRS grid is shared across orbits). Co-Authored-By: Claude Opus 4.8 * docs(s1-rtc): note the dual-orbit time axis spans both orbits in the datacube A dual-orbit cube merges two per-orbit sub-cubes (disjoint time axes) onto a shared grid. Modelling orbit as a dimension would imply a sparse, mostly-empty orbit×time grid; instead orbit stays an attribute of each acquisition, conveyed via the per-orbit assets. Add a `description` on the time dimension so the merged axis isn't misread. Single-orbit cubes get no note. Co-Authored-By: Claude Opus 4.8 * fix(s1-rtc): distinct per-acquisition titles; clearer asset titles - per-acquisition items inherited the cube's "… — tile {id}" title, so every scene in the acquisitions collection was titled identically. Give each a title with its datetime + orbit so siblings are distinguishable. - simplify the border-mask asset title ("Valid-data mask") and make the zarr-store title say "Sentinel-1" for naming consistency. Co-Authored-By: Claude Opus 4.8 --------- Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/stac/s1_rtc.py | 453 ++++++++++++++++++++++---- tests/test_s1_stac.py | 144 ++++++-- tests/test_s1_stac_per_acquisition.py | 286 ++++++++++++++++ 3 files changed, 806 insertions(+), 77 deletions(-) create mode 100644 tests/test_s1_stac_per_acquisition.py diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index eba7da6a..6e50d33e 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -4,7 +4,7 @@ import datetime as dt from pathlib import Path -from typing import cast +from typing import NamedTuple, cast import numpy as np import pyproj @@ -15,8 +15,22 @@ SAT_EXT = "https://stac-extensions.github.io/sat/v1.0.0/schema.json" PROJ_EXT = "https://stac-extensions.github.io/projection/v2.0.0/schema.json" RENDER_EXT = "https://stac-extensions.github.io/render/v1.0.0/schema.json" +DATACUBE_EXT = "https://stac-extensions.github.io/datacube/v2.2.0/schema.json" +TIMESTAMPS_EXT = "https://stac-extensions.github.io/timestamps/v1.1.0/schema.json" + +ZARR_MEDIA_TYPE = "application/vnd.zarr; version=3" _ORBIT_PREFERENCE = ("ascending", "descending") +# Short suffix for orbit-keyed asset names (gamma0-rtc-backscatter-asc / -desc). +_ORBIT_SHORT = {"ascending": "asc", "descending": "desc"} + +# γ⁰ RTC backscatter is float32 with a NaN fill at every resolution level; the arrays carry no attrs +# in the store so these product invariants are hardcoded (see the store hierarchy in data_api/s1_rtc). +GAMMA0_DTYPE = "float32" +GAMMA0_NODATA = "nan" +GAMMA0_UNIT = "gamma0 (linear power)" +BORDER_MASK_DTYPE = "uint8" +GSD = 10 def _rgb_render(orbit: str) -> dict[str, object]: @@ -32,15 +46,28 @@ def _rgb_render(orbit: str) -> dict[str, object]: return { "title": "VV, VH, VV/VH composite", "expression": f"{vv};{vh};({vv})/({vh})", - "rescale": [[0.0, 0.1]], + # Linear gamma0 stretch tuned for the S1 RTC product (the 0.0,0.1 default was too bright). + "rescale": [[0.0, 0.2]], "bidx": [1], "tilesize": 256, } -def _utm_to_wgs84( - proj_code: str, utm_bbox: list[float] -) -> tuple[float, float, float, float]: +def _gamma0_bands() -> list[dict[str, object]]: + """STAC 1.1 band objects for the two polarisations carried by a γ⁰ RTC asset.""" + return [ + { + "name": pol, + "description": f"γ⁰ RTC backscatter, {pol.upper()} polarization", + "data_type": GAMMA0_DTYPE, + "nodata": GAMMA0_NODATA, + "unit": GAMMA0_UNIT, + } + for pol in ("vv", "vh") + ] + + +def _utm_to_wgs84(proj_code: str, utm_bbox: list[float]) -> tuple[float, float, float, float]: """Convert UTM [xmin, ymin, xmax, ymax] to WGS84 (west, south, east, north).""" xmin, ymin, xmax, ymax = utm_bbox transformer = pyproj.Transformer.from_crs(proj_code, "EPSG:4326", always_xy=True) @@ -50,6 +77,33 @@ def _utm_to_wgs84( return min(lons), min(lats), max(lons), max(lats) +def _bbox_to_geometry(bbox: list[float]) -> dict[str, object]: + """A closed rectangular Polygon for a WGS84 [west, south, east, north] bbox.""" + west, south, east, north = bbox + return { + "type": "Polygon", + "coordinates": [ + [[west, south], [east, south], [east, north], [west, north], [west, south]] + ], + } + + +def _open_root(zarr_store: str) -> zarr.Group: + """Open the cube root, preferring consolidated metadata. + + A cube grown by appending a time-slice to an *existing same-orbit* group can end up without root + consolidated metadata (re-consolidating an append on the S3 store is unreliable). The builder must + not require it — fall back to reading the hierarchy directly, exactly as titiler does. See the + data-model issue on the S1 RTC consolidated-metadata regression. + """ + try: + return zarr.open_consolidated(zarr_store, zarr_format=3) + except ValueError as exc: + if "consolidated metadata" not in str(exc).lower(): + raise + return zarr.open_group(zarr_store, mode="r", zarr_format=3) + + def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: """Build a STAC item from a consolidated S1 GRD RTC Zarr V3 store. @@ -74,38 +128,40 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: # href). Revert this prefix to "s1-grd-rtc-" when titiler-eopf#108 lands. tile_id = Path(zarr_store).name.removeprefix("s1-rtc-").removesuffix(".zarr") - # A cube grown by appending a time-slice to an *existing same-orbit* group can end up without - # root consolidated metadata (re-consolidating an append on the S3 store is unreliable). The - # builder must not require it — fall back to reading the hierarchy directly, exactly as titiler - # does. See the data-model issue on the S1 RTC consolidated-metadata regression. - try: - root = zarr.open_consolidated(zarr_store, zarr_format=3) - except ValueError as exc: - if "consolidated metadata" not in str(exc).lower(): - raise - root = zarr.open_group(zarr_store, mode="r", zarr_format=3) + root = _open_root(zarr_store) all_times_ns: list[int] = [] wgs84_bboxes: list[tuple[float, float, float, float]] = [] - preferred_orbit: str | None = None + # Per present orbit, in preference order: the metadata needed for assets, projection and datacube. + present: list[dict[str, object]] = [] for orbit_dir in _ORBIT_PREFERENCE: if orbit_dir not in root: continue - og = cast(zarr.Group, root[orbit_dir]) + og = cast("zarr.Group", root[orbit_dir]) attrs = dict(og.attrs) - proj_code = cast(str, attrs["proj:code"]) - utm_bbox = cast(list[float], attrs["spatial:bbox"]) + proj_code = cast("str", attrs["proj:code"]) + utm_bbox = cast("list[float]", attrs["spatial:bbox"]) - r10m = cast(zarr.Group, og["r10m"]) - times = np.array(cast(zarr.Array, r10m["time"])).tolist() + r10m = cast("zarr.Group", og["r10m"]) + times = np.array(cast("zarr.Array", r10m["time"])).tolist() if not times: continue + # proj:shape / proj:transform live on the r10m group attrs in real stores; read best-effort so + # minimal/legacy stores without them still build (just without those projection refinements). + r10m_attrs = dict(r10m.attrs) all_times_ns.extend(times) wgs84_bboxes.append(_utm_to_wgs84(proj_code, utm_bbox)) - if preferred_orbit is None: - preferred_orbit = orbit_dir + present.append( + { + "orbit": orbit_dir, + "proj_code": proj_code, + "utm_bbox": utm_bbox, + "shape": r10m_attrs.get("spatial:shape"), + "transform": r10m_attrs.get("spatial:transform"), + } + ) if not all_times_ns: raise ValueError(f"No acquisitions found in Zarr store: {zarr_store}") @@ -121,39 +177,113 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: north = max(b[3] for b in wgs84_bboxes) wgs84_bbox = [west, south, east, north] - geometry = { - "type": "Polygon", - "coordinates": [ - [[west, south], [east, south], [east, north], [west, north], [west, south]] - ], + geometry = _bbox_to_geometry(wgs84_bbox) + + # The preferred orbit (ascending if present) drives the single-valued projection fields and the + # default render/preview; every present orbit gets its own first-class asset below. + preferred = present[0] + preferred_orbit = cast("str", preferred["orbit"]) + preferred_proj_code = cast("str", preferred["proj_code"]) + preferred_bbox = cast("list[float]", preferred["utm_bbox"]) + + properties: dict[str, object] = { + "start_datetime": start_dt.isoformat(), + "end_datetime": end_dt.isoformat(), + "title": f"Sentinel-1 GRD RTC γ⁰ — tile {tile_id}", + "description": ( + "Radiometric-terrain-corrected (RTC) γ⁰ backscatter datacube from Sentinel-1 GRD, " + "reprojected onto the Sentinel-2 MGRS/UTM grid." + ), + # `updated` (timestamps extension) tracks this metadata build. `created` is intentionally + # omitted: it means the item's creation instant, which the store does not record — using an + # acquisition time would misuse the field, and a build-time value would churn on every append. + "updated": dt.datetime.now(tz=dt.UTC).isoformat(), + # Identity invariants (constant across the cube; platform is per-acquisition so omitted here — + # a cube can mix S1A and S1C). + "constellation": "sentinel-1", + "instruments": ["c-sar"], + "gsd": GSD, + # SAR extension + "sar:instrument_mode": "IW", + "sar:frequency_band": "C", + "sar:center_frequency": 5.405, + "sar:polarizations": ["VV", "VH"], + "sar:product_type": "GRD", + # Projection extension + "proj:code": preferred_proj_code, + "proj:bbox": preferred_bbox, + # Render extension: dual-pol RGB composite for previews/tiles (defaults to the preferred orbit) + "renders": {"rgb": _rgb_render(preferred_orbit)}, } + if preferred["shape"] is not None: + properties["proj:shape"] = preferred["shape"] + if preferred["transform"] is not None: + properties["proj:transform"] = preferred["transform"] - # preferred_orbit is guaranteed non-None here (ValueError raised above if no acquisitions) - assert preferred_orbit is not None - preferred_proj_code = cast(str, dict(cast(zarr.Group, root[preferred_orbit]).attrs)["proj:code"]) + stac_extensions = [SAR_EXT, PROJ_EXT, RENDER_EXT, DATACUBE_EXT, TIMESTAMPS_EXT] + + # sat:orbit_state is single-valued, so it's only meaningful when the cube holds a single orbit. A + # dual-orbit cube would mislabel half its slices — omit it there (per-acquisition items, which are + # single-orbit, carry the real value). Only declare the SAT extension when the field is set. + if len(present) == 1: + properties["sat:orbit_state"] = preferred_orbit + stac_extensions.append(SAT_EXT) + + # Datacube extension. The time axis is irregularly sampled, so it lists its discrete `values` (the + # acquisition instants across all orbit groups, sorted) — their count is the number of time steps, + # and the list stays modest (bounded by the tile's acquisitions). The regular x/y axes instead carry + # extent + step (their element count is derivable, and the exact pixel count is in proj:shape); + # enumerating their ~10⁴ coordinates would not scale. + epsg = pyproj.CRS.from_user_input(preferred_proj_code).to_epsg() + xmin, ymin, xmax, ymax = preferred_bbox + time_values = [ + dt.datetime.fromtimestamp(t / 1e9, tz=dt.UTC).isoformat() for t in sorted(set(all_times_ns)) + ] + time_dim: dict[str, object] = { + "type": "temporal", + "extent": [start_dt.isoformat(), end_dt.isoformat()], + "values": time_values, + } + if len(present) > 1: + # The cube merges two per-orbit sub-cubes (disjoint time axes) onto a shared grid. Orbit is an + # attribute of each acquisition, not an independent axis, so it is conveyed via the per-orbit + # assets rather than a (sparse) orbit dimension — note that here to avoid misreading the axis. + time_dim["description"] = ( + "Acquisition instants across both orbit directions (union); each step belongs to a single " + "orbit — the ascending/descending groups are exposed as separate assets and as items in " + "the per-acquisition collection." + ) + x_dim: dict[str, object] = { + "type": "spatial", + "axis": "x", + "extent": [xmin, xmax], + "reference_system": epsg, + } + y_dim: dict[str, object] = { + "type": "spatial", + "axis": "y", + "extent": [ymin, ymax], + "reference_system": epsg, + } + if preferred["transform"] is not None: + transform = cast("list[float]", preferred["transform"]) + x_dim["step"] = transform[0] + y_dim["step"] = transform[4] + properties["cube:dimensions"] = {"time": time_dim, "x": x_dim, "y": y_dim} + # `variable_type` is the datacube field name (not `type`); the border mask is auxiliary, not data. + properties["cube:variables"] = { + "vv": {"dimensions": ["time", "y", "x"], "variable_type": "data", "unit": GAMMA0_UNIT}, + "vh": {"dimensions": ["time", "y", "x"], "variable_type": "data", "unit": GAMMA0_UNIT}, + "border_mask": {"dimensions": ["time", "y", "x"], "variable_type": "auxiliary"}, + } item = pystac.Item( id=f"s1-rtc-{tile_id}", geometry=geometry, bbox=wgs84_bbox, datetime=None, - properties={ - "start_datetime": start_dt.isoformat(), - "end_datetime": end_dt.isoformat(), - # SAR extension - "sar:instrument_mode": "IW", - "sar:frequency_band": "C", - "sar:center_frequency": 5.405, - "sar:polarizations": ["VV", "VH"], - "sar:product_type": "GRD", - # SAT extension - "sat:orbit_state": preferred_orbit, - # Projection extension - "proj:code": preferred_proj_code, - # Render extension: dual-pol RGB composite for previews/tiles - "renders": {"rgb": _rgb_render(preferred_orbit)}, - }, - stac_extensions=[SAR_EXT, SAT_EXT, PROJ_EXT, RENDER_EXT], + properties=properties, + stac_extensions=stac_extensions, collection=collection_id, ) @@ -162,20 +292,233 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: "zarr-store", pystac.Asset( href=store_str, - media_type="application/vnd.zarr; version=3", + media_type=ZARR_MEDIA_TYPE, roles=["data"], - title="S1 GRD RTC Zarr store", + title="Sentinel-1 GRD RTC Zarr store", ), ) - for pol in ("vv", "vh"): + + # One γ⁰ asset per present orbit group (fixes the duplicate-href vv/vh ambiguity and the missing + # descending asset): VV/VH are addressable as named `bands`, not indistinguishable duplicate assets. + # A separate border-mask asset exposes the valid-data mask variable in the same group. + for info in present: + orbit = cast("str", info["orbit"]) + short = _ORBIT_SHORT[orbit] + group_href = f"{store_str}/{orbit}" item.add_asset( - pol, + f"gamma0-rtc-backscatter-{short}", pystac.Asset( - href=f"{store_str}/{preferred_orbit}", - media_type="application/vnd.zarr; version=3", + href=group_href, + media_type=ZARR_MEDIA_TYPE, roles=["data"], - title=f"{pol.upper()} polarisation", + title=f"γ⁰ RTC backscatter ({orbit})", + extra_fields={ + "bands": _gamma0_bands(), + "data_type": GAMMA0_DTYPE, + "nodata": GAMMA0_NODATA, + "unit": GAMMA0_UNIT, + "gsd": GSD, + }, + ), + ) + item.add_asset( + f"border-mask-{short}", + pystac.Asset( + href=group_href, + media_type=ZARR_MEDIA_TYPE, + roles=["data"], + title=f"Valid-data mask ({orbit})", + extra_fields={ + "bands": [ + { + "name": "border_mask", + "description": "Valid-data mask (0 = border/no-data, non-zero = valid)", + "data_type": BORDER_MASK_DTYPE, + "nodata": 0, + } + ], + "gsd": GSD, + }, ), ) return item + + +# ============================================================================ +# Per-acquisition item construction (one queryable item per cube `time` slice) +# ============================================================================ + +# Default the cube preview to the most recent acquisition covering most of the tile, so a browser shows +# fresh near-full data rather than the oldest slice. +COVERAGE_THRESHOLD = 0.80 + + +class Slice(NamedTuple): + """One cube time slice: its orbit group, acquisition instant, and tile coverage fraction (0..1).""" + + orbit: str + dt: dt.datetime + coverage: float + + +def pick_slice(slices: list[Slice]) -> Slice | None: + """Choose the slice the cube preview should default to. + + The most recent acquisition with coverage strictly above ``COVERAGE_THRESHOLD``; if none clears it, + the highest-coverage slice (ties broken by most recent). Spans both orbit groups. Returns ``None`` + for an empty cube. + """ + if not slices: + return None + good = [s for s in slices if s.coverage > COVERAGE_THRESHOLD] + if good: + return max(good, key=lambda s: s.dt) + return max(slices, key=lambda s: (s.coverage, s.dt)) + + +def slice_coverages(zarr_store: str) -> list[Slice]: + """Per-slice tile coverage from the cube, across both orbit groups. + + Reads ``border_mask`` at the cheap ``r720m`` overview only (~150x150). Coverage is the fraction of + **valid** pixels; the S1Tiling border mask is stored with ``fill_value=0`` for the border, so valid + = non-zero. ``time`` is raw int64 ns (as :func:`build_s1_rtc_stac_item` reads it) -> UTC datetime. + """ + root = _open_root(zarr_store) + out: list[Slice] = [] + for orbit in _ORBIT_PREFERENCE: + if orbit not in root: + continue + level = cast("zarr.Group", cast("zarr.Group", root[orbit])["r720m"]) + mask = np.asarray(cast("zarr.Array", level["border_mask"])) # (time, y, x), uint8 + times_ns = np.asarray(cast("zarr.Array", level["time"])).tolist() # int64 ns since epoch + for i, t_ns in enumerate(times_ns): + sl = mask[i] + coverage = float(np.count_nonzero(sl) / sl.size) + out.append(Slice(orbit, dt.datetime.fromtimestamp(t_ns / 1e9, tz=dt.UTC), coverage)) + return out + + +def acquisition_id(tile_id: str, when: dt.datetime) -> str: + """Per-acquisition item id, e.g. ``s1-rtc-31TCH-20260607t055248``.""" + return f"s1-rtc-{tile_id}-{when.strftime('%Y%m%dt%H%M%S')}" + + +def _normalize_platform(raw: object) -> str | None: + """Map the store's short platform code (e.g. ``s1a``) to the STAC convention (``sentinel-1a``). + + Mirrors the Sentinel-2 reference (``sentinel-2a``). Unknown values are returned lower-cased. + """ + s = str(raw).strip().lower() + if not s: + return None + if len(s) == 3 and s.startswith("s1"): + return f"sentinel-1{s[2]}" + return s + + +def build_s1_rtc_per_acquisition_items( + zarr_store: str, *, orbit: str, collection_id: str +) -> list[pystac.Item]: + """Build one queryable STAC item per cube ``time`` slice of a single orbit group. + + Each item is a single-``datetime`` view into the shared cube (no data duplication): it keeps the + cube's geometry/SAR/projection metadata and the orbit's γ⁰ asset, drops the temporal range and the + datacube structure (a single acquisition is not a cube), and is reoriented to ``orbit``. The item + is deployment-agnostic — it carries the render config + datetime, and the registration layer derives + the TiTiler links (which point at the cube endpoint with ``sel=time={datetime}``) from it. + + Parameters + ---------- + zarr_store: + Local path or ``s3://`` URI to the per-tile cube Zarr store. + orbit: + Orbit group to emit items for (``"ascending"`` or ``"descending"``). + collection_id: + Target (per-acquisition) STAC collection ID. + + Raises + ------ + ValueError + If the store has no acquisitions, or ``orbit`` is not present in the store. + """ + if orbit not in _ORBIT_PREFERENCE: + raise ValueError(f"orbit must be one of {_ORBIT_PREFERENCE}, got {orbit!r}") + + tile_id = Path(zarr_store).name.removeprefix("s1-rtc-").removesuffix(".zarr") + + root = _open_root(zarr_store) + if orbit not in root: + raise ValueError(f"Orbit group {orbit!r} not found in Zarr store: {zarr_store}") + r10m = cast("zarr.Group", cast("zarr.Group", root[orbit])["r10m"]) + times_ns = np.array(cast("zarr.Array", r10m["time"])).tolist() + platforms = np.array(cast("zarr.Array", r10m["platform"])).tolist() + if not times_ns: + raise ValueError(f"No acquisitions found for orbit {orbit!r} in: {zarr_store}") + + base = build_s1_rtc_stac_item(zarr_store, collection_id) + base_dict = base.to_dict(include_self_link=False) + + # Assets to drop from each per-acquisition clone: the *other* orbit's groups (a per-acq item + # represents one orbit). The datacube structure is dropped too — a single acquisition is not a cube. + other_assets = { + key + for o in _ORBIT_PREFERENCE + if o != orbit + for key in (f"gamma0-rtc-backscatter-{_ORBIT_SHORT[o]}", f"border-mask-{_ORBIT_SHORT[o]}") + } + + # A per-acquisition item covers only its run orbit's footprint — not the cube's union of both + # orbits' extents (which the base item carries). Recompute bbox/geometry/proj:bbox from this orbit. + # proj:code/shape/transform describe the shared MGRS grid (identical across orbits), so the values + # inherited from the base (preferred orbit) are correct and are intentionally not recomputed here. + og_attrs = dict(cast("zarr.Group", root[orbit]).attrs) + orbit_utm_bbox = cast("list[float]", og_attrs["spatial:bbox"]) + orbit_wgs84_bbox = list(_utm_to_wgs84(cast("str", og_attrs["proj:code"]), orbit_utm_bbox)) + orbit_geometry = _bbox_to_geometry(orbit_wgs84_bbox) + + items: list[pystac.Item] = [] + for t_ns, platform in zip(times_ns, platforms, strict=True): + when = dt.datetime.fromtimestamp(t_ns / 1e9, tz=dt.UTC) + item_dict = {**base_dict} + item_dict["id"] = acquisition_id(tile_id, when) + item_dict["bbox"] = orbit_wgs84_bbox + item_dict["geometry"] = orbit_geometry + + props = { + k: v + for k, v in base_dict["properties"].items() + if k not in ("start_datetime", "end_datetime", "cube:dimensions", "cube:variables") + } + props["datetime"] = when.isoformat() + props["sat:orbit_state"] = orbit + props["proj:bbox"] = orbit_utm_bbox + # Per-acquisition title carries the datetime + orbit so sibling scenes are distinguishable + # (the inherited cube title "… — tile {id}" is identical across all acquisitions). + props["title"] = ( + f"Sentinel-1 GRD RTC γ⁰ — tile {tile_id}, " + f"{when.strftime('%Y-%m-%dT%H:%M:%SZ')} ({orbit})" + ) + props["description"] = ( + "Radiometric-terrain-corrected (RTC) γ⁰ backscatter from a single Sentinel-1 GRD " + "acquisition, reprojected onto the Sentinel-2 MGRS/UTM grid." + ) + normalized = _normalize_platform(platform) + if normalized: + props["platform"] = normalized + props["renders"] = {"rgb": _rgb_render(orbit)} + item_dict["properties"] = props + + # Drop the datacube ext (a single acquisition is not a cube). Ensure the SAT ext is declared: + # a per-acq item always sets sat:orbit_state, but a dual-orbit cube base omits both. + exts = [e for e in base_dict.get("stac_extensions", []) if e != DATACUBE_EXT] + if SAT_EXT not in exts: + exts.append(SAT_EXT) + item_dict["stac_extensions"] = exts + item_dict["assets"] = { + k: v for k, v in base_dict["assets"].items() if k not in other_assets + } + item_dict["links"] = [] + + items.append(pystac.Item.from_dict(item_dict)) + return items diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index cb2b9226..732d880f 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -52,16 +52,29 @@ def _make_s1_store( # render path resolves; revert to "s1-grd-rtc-" when titiler-eopf#108 lands. store_path = tmp_path / f"s1-rtc-{tile_id}.zarr" root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + ny = nx = 4 # tiny spatial grid: enough for the builder's metadata reads for orbit_dir, acquisitions in orbits.items(): og = root.create_group(orbit_dir) og.attrs.update({"proj:code": crs, "spatial:bbox": utm_bbox}) r10m = og.create_group("r10m") + # proj:shape / proj:transform live on the r10m group attrs in real stores. + r10m.attrs.update( + { + "spatial:shape": [ny, nx], + "spatial:transform": [10.0, 0.0, utm_bbox[0], 0.0, -10.0, utm_bbox[3]], + } + ) times = np.array([t for t, _ in acquisitions], dtype="int64") platforms = np.array([p for _, p in acquisitions], dtype=" None: # Union must be wider than either individual bbox west, _south, east, _north = item.bbox # type: ignore[misc] - assert west < 1.0 # left edge from descending - assert east > 2.5 # right edge from ascending (shifted ~1° further east) + assert west < 1.0 # left edge from descending + assert east > 2.5 # right edge from ascending (shifted ~1° further east) -def test_ascending_preferred_for_assets(tmp_path: Path) -> None: - """When both orbits present, vv/vh assets must point to the ascending group.""" - store_path = tmp_path / "s1-rtc-31TCH.zarr" - root = zarr.open_group(str(store_path), mode="w", zarr_format=3) - for orbit_dir in ("descending", "ascending"): - og = root.create_group(orbit_dir) - og.attrs.update({"proj:code": CRS, "spatial:bbox": UTM_BBOX}) - r10m = og.create_group("r10m") - t_arr = r10m.create_array("time", shape=(1,), dtype="int64", chunks=(512,)) - t_arr[:] = [T1_NS] - p_arr = r10m.create_array("platform", shape=(1,), dtype=" None: + """A dual-orbit cube must expose a γ⁰ asset per orbit group (bug #2), each pointing at its group.""" + store = _make_s1_store( + tmp_path, {"descending": [(T1_NS, "S1A")], "ascending": [(T1_NS, "S1A")]} + ) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") - item = build_s1_rtc_stac_item(str(store_path), "sentinel-1-grd-rtc-staging") - assert "ascending" in item.assets["vv"].href - assert "ascending" in item.assets["vh"].href + asc = item.assets["gamma0-rtc-backscatter-asc"] + desc = item.assets["gamma0-rtc-backscatter-desc"] + assert asc.href.endswith("/ascending") + assert desc.href.endswith("/descending") + # The dual-pol href ambiguity (bug #1) is gone: VV/VH are named bands, not duplicate assets. + assert [b["name"] for b in asc.extra_fields["bands"]] == ["vv", "vh"] + assert "border-mask-asc" in item.assets + assert "border-mask-desc" in item.assets def test_empty_store_raises(tmp_path: Path) -> None: @@ -188,14 +199,103 @@ def test_empty_store_raises(tmp_path: Path) -> None: def test_asset_hrefs(tmp_path: Path) -> None: - """zarr-store href = store URI; vv/vh hrefs = {store}/{orbit} (orbit group root, per geozarr spec).""" + """zarr-store href = store URI; the γ⁰ asset href = {store}/{orbit} (orbit group, per geozarr spec).""" store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") store_str = str(store) assert item.assets["zarr-store"].href == store_str - assert item.assets["vv"].href == f"{store_str}/descending" - assert item.assets["vh"].href == f"{store_str}/descending" + # Single-orbit cube → only the descending γ⁰/mask assets exist (no ascending). + assert item.assets["gamma0-rtc-backscatter-desc"].href == f"{store_str}/descending" + assert item.assets["border-mask-desc"].href == f"{store_str}/descending" + assert "gamma0-rtc-backscatter-asc" not in item.assets + + +def test_gamma0_asset_band_and_invariant_metadata(tmp_path: Path) -> None: + """The γ⁰ asset carries VV/VH bands + data_type/nodata/unit/gsd invariants.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + asset = item.assets["gamma0-rtc-backscatter-desc"] + assert asset.extra_fields["data_type"] == "float32" + assert asset.extra_fields["nodata"] == "nan" + assert asset.extra_fields["gsd"] == 10 + bands = asset.extra_fields["bands"] + assert {b["name"] for b in bands} == {"vv", "vh"} + assert all(b["data_type"] == "float32" for b in bands) + + +def test_identity_and_projection_fields(tmp_path: Path) -> None: + """Item-level identity invariants + projection detail are present (S2-parity gaps).""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + props = item.properties + assert props["constellation"] == "sentinel-1" + assert props["instruments"] == ["c-sar"] + assert props["gsd"] == 10 + assert "platform" not in props # per-acquisition; a cube can mix S1A/S1C + assert props["proj:bbox"] == UTM_BBOX + assert props["proj:shape"] == [4, 4] + assert props["proj:transform"][0] == 10.0 + + +def test_datacube_extension(tmp_path: Path) -> None: + """Cube items carry the datacube extension. The irregular time axis lists its discrete `values` + (count = number of acquisitions); the regular x/y axes carry extent + step only (no per-pixel + enumeration — the exact pixel count is in proj:shape).""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A"), (T2_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + assert "https://stac-extensions.github.io/datacube/v2.2.0/schema.json" in item.stac_extensions + dims = item.properties["cube:dimensions"] + assert dims["time"]["type"] == "temporal" + assert len(dims["time"]["extent"]) == 2 + assert dims["time"]["values"] == [ # two distinct acquisitions -> two time steps, sorted + dt.datetime.fromtimestamp(T1_NS / 1e9, tz=dt.UTC).isoformat(), + dt.datetime.fromtimestamp(T2_NS / 1e9, tz=dt.UTC).isoformat(), + ] + assert dims["x"]["reference_system"] == 32631 + assert dims["x"]["step"] == 10.0 + assert dims["y"]["step"] == -10.0 + assert "values" not in dims["x"] # regular axis: extent + step, not ~10^4 coordinates + assert item.properties["proj:shape"] == [4, 4] # exact x/y element count + variables = item.properties["cube:variables"] + assert set(variables) == {"vv", "vh", "border_mask"} + assert variables["vv"]["dimensions"] == ["time", "y", "x"] + # datacube field is `variable_type` (not `type`); the border mask is auxiliary, not data. + assert variables["vv"]["variable_type"] == "data" + assert variables["border_mask"]["variable_type"] == "auxiliary" + assert "type" not in variables["vv"] + + +def test_orbit_state_single_vs_dual(tmp_path: Path) -> None: + """sat:orbit_state (single-valued) is set only for a single-orbit cube; a dual-orbit cube omits it + (and the SAT extension) rather than mislabel half its slices.""" + single = _make_s1_store(tmp_path / "a", {"descending": [(T1_NS, "S1A")]}) + item1 = build_s1_rtc_stac_item(str(single), "sentinel-1-grd-rtc-staging") + assert item1.properties["sat:orbit_state"] == "descending" + assert "https://stac-extensions.github.io/sat/v1.0.0/schema.json" in item1.stac_extensions + assert "description" not in item1.properties["cube:dimensions"]["time"] # single orbit: no note + + dual = _make_s1_store( + tmp_path / "b", {"ascending": [(T1_NS, "S1A")], "descending": [(T1_NS, "S1A")]} + ) + item2 = build_s1_rtc_stac_item(str(dual), "sentinel-1-grd-rtc-staging") + assert "sat:orbit_state" not in item2.properties + assert "https://stac-extensions.github.io/sat/v1.0.0/schema.json" not in item2.stac_extensions + # dual-orbit cube notes that the merged time axis spans both orbits (orbit is an asset-level split) + assert "orbit" in item2.properties["cube:dimensions"]["time"]["description"].lower() + + +def test_timestamps_updated_only(tmp_path: Path) -> None: + """`updated` (build time) is set; `created` is omitted — the store records no item-creation time, + so an acquisition time would misuse the field and a build-time value would churn on every append.""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A"), (T2_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + + assert "created" not in item.properties + assert dt.datetime.fromisoformat(item.properties["updated"]) def test_sar_extension_fields(tmp_path: Path) -> None: @@ -224,7 +324,7 @@ def test_render_extension_rgb_composite(tmp_path: Path) -> None: rgb = item.properties["renders"]["rgb"] assert rgb["expression"] == "/descending:vv;/descending:vh;(/descending:vv)/(/descending:vh)" - assert rgb["rescale"] == [[0.0, 0.1]] + assert rgb["rescale"] == [[0.0, 0.2]] assert rgb["bidx"] == [1] assert rgb["tilesize"] == 256 diff --git a/tests/test_s1_stac_per_acquisition.py b/tests/test_s1_stac_per_acquisition.py new file mode 100644 index 00000000..8411e9eb --- /dev/null +++ b/tests/test_s1_stac_per_acquisition.py @@ -0,0 +1,286 @@ +"""Tests for the per-acquisition + coverage construction moved into eopf_geozarr.stac.s1_rtc. + +Ported from the data-pipeline construction tests (``test_pick_slice``, ``test_slice_coverage`` and the +construction subset of ``test_register_per_acquisition``). The TiTiler-link / thumbnail / alternate-asset +behaviour stays in data-pipeline — those are registration concerns, not item construction. +""" + +from __future__ import annotations + +import datetime as dt +import json +from typing import TYPE_CHECKING + +import numpy as np +import pytest +import zarr + +from eopf_geozarr.stac.s1_rtc import ( + DATACUBE_EXT, + Slice, + acquisition_id, + build_s1_rtc_per_acquisition_items, + pick_slice, + slice_coverages, +) + +if TYPE_CHECKING: + from pathlib import Path + +CRS = "EPSG:32631" +UTM_BBOX = [300000.0, 4900000.0, 400000.0, 5000000.0] + +# Two acquisitions in PHYSICAL (append) order — deliberately NOT chronological, to prove each item's +# datetime follows its own slice, not its list position. +T_LATER = int(dt.datetime(2026, 6, 7, 5, 52, 48, tzinfo=dt.UTC).timestamp() * 1e9) +T_EARLY = int(dt.datetime(2026, 6, 5, 6, 9, 7, tzinfo=dt.UTC).timestamp() * 1e9) + + +# ============================================================================= +# pick_slice — pure preview-slice selector +# ============================================================================= + + +def _s(orbit: str, day: int, coverage: float) -> Slice: + return Slice(orbit=orbit, dt=dt.datetime(2026, 6, day, 6, 0, tzinfo=dt.UTC), coverage=coverage) + + +def test_most_recent_above_threshold_wins_even_if_older_has_more_coverage() -> None: + chosen = pick_slice( + [_s("ascending", 4, 0.99), _s("descending", 7, 0.85), _s("ascending", 6, 0.90)] + ) + assert chosen is not None + assert (chosen.orbit, chosen.dt.day, chosen.coverage) == ("descending", 7, 0.85) + + +def test_falls_back_to_max_coverage_when_none_above_threshold() -> None: + chosen = pick_slice( + [_s("ascending", 7, 0.40), _s("descending", 5, 0.75), _s("ascending", 4, 0.50)] + ) + assert chosen is not None + assert (chosen.orbit, chosen.dt.day, chosen.coverage) == ("descending", 5, 0.75) + + +def test_exact_80_percent_is_not_above_threshold() -> None: + # 0.80 is NOT > 0.80 -> excluded from the "good" set; the 0.81 older slice wins instead. + chosen = pick_slice([_s("ascending", 7, 0.80), _s("descending", 4, 0.81)]) + assert chosen is not None + assert (chosen.dt.day, chosen.coverage) == (4, 0.81) + + +def test_max_coverage_tie_broken_by_most_recent() -> None: + chosen = pick_slice( + [_s("ascending", 4, 0.60), _s("descending", 8, 0.60), _s("ascending", 6, 0.60)] + ) + assert chosen is not None + assert chosen.dt.day == 8 + + +def test_single_slice_returned() -> None: + only = _s("ascending", 5, 0.10) + assert pick_slice([only]) == only + + +def test_empty_returns_none() -> None: + assert pick_slice([]) is None + + +# ============================================================================= +# slice_coverages — per-slice tile coverage from the cube's r720m border_mask +# ============================================================================= + + +def _ns(day: int) -> int: + return int(dt.datetime(2026, 6, day, 6, 0, tzinfo=dt.UTC).timestamp() * 1e9) + + +def _write_r720m(root: zarr.Group, orbit: str, masks: list[np.ndarray], days: list[int]) -> None: + """Write {orbit}/r720m with a border_mask (time,y,x) + time (int64 ns), like the real cube.""" + lvl = root.create_group(orbit).create_group("r720m") + n = len(masks) + y, x = masks[0].shape + bm = lvl.create_array("border_mask", shape=(n, y, x), dtype="uint8", fill_value=0) + bm[:] = np.stack(masks).astype("uint8") + t = lvl.create_array("time", shape=(n,), dtype="int64") + t[:] = np.array([_ns(d) for d in days], dtype="int64") + + +def _make_coverage_cube(tmp_path: Path) -> str: + store = str(tmp_path / "s1-rtc-31TCH.zarr") + root = zarr.open_group(store, mode="w", zarr_format=3) + full = np.ones((4, 4)) # coverage 1.0 + half = np.zeros((4, 4)) + half[:2, :] = 1 # coverage 0.5 + quarter = np.zeros((4, 4)) + quarter[0, :] = 1 # coverage 0.25 + _write_r720m(root, "ascending", [full, half], [4, 6]) + _write_r720m(root, "descending", [quarter], [5]) + return store + + +def test_slice_coverages_reads_both_orbits_with_correct_fractions(tmp_path: Path) -> None: + by_key = { + (s.orbit, s.dt.day): s.coverage for s in slice_coverages(_make_coverage_cube(tmp_path)) + } + assert by_key == { + ("ascending", 4): 1.0, # full mask -> polarity check (non-zero = valid) + ("ascending", 6): 0.5, + ("descending", 5): 0.25, + } + + +def test_slice_coverages_times_are_utc_datetimes(tmp_path: Path) -> None: + s = next(s for s in slice_coverages(_make_coverage_cube(tmp_path)) if s.orbit == "descending") + assert s.dt == dt.datetime(2026, 6, 5, 6, 0, tzinfo=dt.UTC) + + +def test_slice_coverages_skips_missing_orbit(tmp_path: Path) -> None: + store = str(tmp_path / "s1-rtc-asc-only.zarr") + root = zarr.open_group(store, mode="w", zarr_format=3) + _write_r720m(root, "ascending", [np.ones((4, 4))], [7]) + assert {s.orbit for s in slice_coverages(store)} == {"ascending"} + + +# ============================================================================= +# build_s1_rtc_per_acquisition_items — one item per cube time slice +# ============================================================================= + + +def _make_acq_cube( + tmp_path: Path, orbits: dict[str, list[tuple[int, str]]], tile_id: str = "31TCH" +) -> str: + """A cube with r10m time/platform (+ tiny data arrays) per orbit — enough to build per-acq items.""" + store = str(tmp_path / f"s1-rtc-{tile_id}.zarr") + root = zarr.open_group(store, mode="w", zarr_format=3) + ny = nx = 4 + for orbit, acqs in orbits.items(): + og = root.create_group(orbit) + og.attrs.update({"proj:code": CRS, "spatial:bbox": UTM_BBOX}) + r10m = og.create_group("r10m") + r10m.attrs.update( + { + "spatial:shape": [ny, nx], + "spatial:transform": [10.0, 0.0, UTM_BBOX[0], 0.0, -10.0, UTM_BBOX[3]], + } + ) + times = np.array([t for t, _ in acqs], dtype="int64") + platforms = np.array([p for _, p in acqs], dtype=" None: + when = dt.datetime(2026, 6, 7, 5, 52, 48, tzinfo=dt.UTC) + assert acquisition_id("31TCH", when) == "s1-rtc-31TCH-20260607t055248" + + +def test_per_acq_title_carries_datetime_and_orbit(tmp_path: Path) -> None: + """Per-acq titles must distinguish sibling scenes (not inherit the cube's '— tile {id}' title).""" + store = _make_acq_cube(tmp_path, {"descending": [(T_EARLY, "s1a")]}) + item = build_s1_rtc_per_acquisition_items(store, orbit="descending", collection_id="acq")[0] + title = item.properties["title"] + assert "2026-06-05T06:09:07Z" in title + assert "descending" in title + assert title != "Sentinel-1 GRD RTC γ⁰ — tile 31TCH" + + +def test_single_datetime_no_range_and_targets_collection(tmp_path: Path) -> None: + store = _make_acq_cube(tmp_path, {"descending": [(T_EARLY, "S1A")]}) + item = build_s1_rtc_per_acquisition_items( + store, orbit="descending", collection_id="sentinel-1-grd-rtc-acquisitions" + )[0] + props = item.properties + assert item.collection_id == "sentinel-1-grd-rtc-acquisitions" + assert props["datetime"] == "2026-06-05T06:09:07+00:00" + assert "start_datetime" not in props + assert "end_datetime" not in props + assert props["sar:product_type"] == "GRD" + + +def test_reorients_orbit_metadata_and_keeps_only_run_orbit_asset(tmp_path: Path) -> None: + """A descending run carries /descending metadata + only the descending γ⁰/mask assets — never the + cube's preferred /ascending.""" + store = _make_acq_cube( + tmp_path, {"ascending": [(T_EARLY, "S1A")], "descending": [(T_EARLY, "S1A")]} + ) + item = build_s1_rtc_per_acquisition_items(store, orbit="descending", collection_id="acq")[0] + assert item.properties["sat:orbit_state"] == "descending" + # SAT ext must be declared even though the dual-orbit cube base omitted it. + assert "https://stac-extensions.github.io/sat/v1.0.0/schema.json" in item.stac_extensions + assert item.properties["renders"]["rgb"]["expression"].startswith("/descending:vv") + assert "gamma0-rtc-backscatter-desc" in item.assets + assert "gamma0-rtc-backscatter-asc" not in item.assets + assert "border-mask-asc" not in item.assets + assert "ascending" not in json.dumps(item.to_dict(include_self_link=False)) + + +def test_per_slice_platform_and_no_datacube(tmp_path: Path) -> None: + """Each item gets its own slice's platform (normalized to the STAC convention); a single + acquisition is not a datacube, and `created` is not set.""" + store = _make_acq_cube(tmp_path, {"descending": [(T_LATER, "s1a"), (T_EARLY, "s1c")]}) + items = build_s1_rtc_per_acquisition_items(store, orbit="descending", collection_id="acq") + assert [i.properties["platform"] for i in items] == ["sentinel-1a", "sentinel-1c"] + for item in items: + assert "cube:dimensions" not in item.properties + assert DATACUBE_EXT not in item.stac_extensions + assert "created" not in item.properties + + +def test_footprint_is_run_orbit_not_cube_union(tmp_path: Path) -> None: + """A per-acq item must carry ITS orbit's footprint, not the cube's union of both orbits' extents.""" + # ascending shifted east so the cube union bbox is wider than the descending orbit alone. + store = str(tmp_path / "s1-rtc-31TCH.zarr") + root = zarr.open_group(store, mode="w", zarr_format=3) + ny = nx = 4 + desc_bbox = [300000.0, 4900000.0, 400000.0, 5000000.0] + asc_bbox = [400000.0, 4900000.0, 500000.0, 5000000.0] + for orbit, bbox in (("descending", desc_bbox), ("ascending", asc_bbox)): + og = root.create_group(orbit) + og.attrs.update({"proj:code": CRS, "spatial:bbox": bbox}) + r10m = og.create_group("r10m") + r10m.attrs.update( + { + "spatial:shape": [ny, nx], + "spatial:transform": [10.0, 0.0, bbox[0], 0.0, -10.0, bbox[3]], + } + ) + r10m.create_array("time", shape=(1,), dtype="int64", chunks=(512,))[:] = [T_EARLY] + r10m.create_array("platform", shape=(1,), dtype=" None: + store = _make_acq_cube(tmp_path, {"descending": [(T_EARLY, "S1A")]}) + with pytest.raises(ValueError, match="orbit must be one of"): + build_s1_rtc_per_acquisition_items(store, orbit="sideways", collection_id="acq") + + +def test_orbit_absent_from_store_raises(tmp_path: Path) -> None: + store = _make_acq_cube(tmp_path, {"descending": [(T_EARLY, "S1A")]}) + with pytest.raises(ValueError, match="not found"): + build_s1_rtc_per_acquisition_items(store, orbit="ascending", collection_id="acq") + + +def test_datetime_follows_slice_not_physical_position(tmp_path: Path) -> None: + """Items are emitted in physical (append) order, each carrying its OWN datetime — so a + non-monotonic cube still yields a correct item per slice.""" + store = _make_acq_cube(tmp_path, {"descending": [(T_LATER, "S1A"), (T_EARLY, "S1A")]}) + items = build_s1_rtc_per_acquisition_items(store, orbit="descending", collection_id="acq") + assert items[0].properties["datetime"] == "2026-06-07T05:52:48+00:00" + assert items[1].properties["datetime"] == "2026-06-05T06:09:07+00:00" + assert items[0].id == "s1-rtc-31TCH-20260607t055248" From 993b69c2d0b7ec7e61091df14d09449ac02da76e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Fri, 19 Jun 2026 10:33:47 +0100 Subject: [PATCH 12/19] perf(s1-rtc): shard the conditions arrays (gamma_area/LIA) like the vv/vh pyramid (#197) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * perf(s1-rtc): shard the conditions arrays like the vv/vh pyramid A real S1 RTC cube is 3807 objects / 3.5 GB, of which 3604 (94.7%) are the conditions/gamma_area_ arrays: [10980,10980] float32, inner chunk 366², NO sharding_indexed -> ~900 tiny chunk objects each. They are time-invariant yet dominate the object count, which dominated the ingest's S3 upload wall-time (a live pod sat ~34 min in "Uploading store" at 9 millicores). The vv/vh/border_mask display pyramid is already sharded (one shard per time slice over the full spatial extent, inner 366²). Apply that same existing layout to the condition arrays: add shards=(h, w) to the one create_array in ingest_s1tiling_conditions. All condition arrays (gamma_area, lia, incidence_angle) share that write path and the same 2D full-resolution shape, so all collapse from ~900 chunk objects to 1 shard object (cube ~3807 -> ~210). calculate_aligned_chunk_size returns a divisor of the dimension, so (h, w) is a clean multiple of the inner chunk (Zarr v3 shard-divisibility). conditions arrays are NOT in TiTiler's render path (vv/vh/border_mask), so this does not touch the web-render layout; it only makes a client read a condition array in one ranged GET instead of ~900. Values are byte-identical. Tests: +2 (sharding codec present; 9 inner chunks -> 1 on-disk shard object + byte-identical roundtrip). 57 passed. Spec: claude-docs/specs/s1_gamma_area_sharding.md. Cross-repo Task T5 of data-pipeline/claude-docs/plans/s1_ingest_upload_perf.md. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP * docs(s1-rtc): record real-S3 sharding benchmark in the T5 spec Validated on the live OVH bucket (laptop->DE): object collapse 100->1 (prod ~900->1), PUT 1.7x faster even with batched concurrency on, divisibility valid at the production 10980² (aligned 366 divides 10980), reads byte-identical. Honest caveat recorded: a full-array sequential read is NOT faster sharded (same bytes, one un-parallelizable object) — the win is object-count (upload + listing) and windowed/partial cloud reads, not full-read throughput. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP * refactor(s1-rtc): trim the conditions-sharding comment to match vv/vh style Comment-only. The surrounding vv/vh sharding has no inline explainer; the long cloud-access rationale now lives in claude-docs/specs/s1_gamma_area_sharding.md. Keep only the non-obvious bits: why one shard, and the Zarr v3 shard-divisibility invariant. No behavior change (20 condition/shard tests green). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP * chore(s1-rtc): drop the claude-docs spec from this PR The data-model repo has no claude-docs/specs convention; the spec was noise for this PR's reviewers. The problem statement, real-S3 benchmark and migration note live in the data-pipeline plan + tracking issue EOPF-Explorer/data-pipeline#288 and PR #197's description. Also drop the now-dangling spec path from the code comment (the rationale stays inline). No behavior change (20 condition/shard tests green). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP --------- Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/conversion/s1_ingest.py | 14 +++++-- tests/test_s1_rtc_ingest.py | 51 ++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 4 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 751b0913..62d723e0 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -985,14 +985,20 @@ def ingest_s1tiling_conditions( conditions[array_name][:, :] = data log.info("Overwrote condition array", array_name=array_name) else: + # Shard like the vv/vh pyramid: one shard over the full (y, x) extent so a 10980² + # condition array is a single object, not ~900 tiny 366²-chunk objects. + # calculate_aligned_chunk_size returns a divisor of the dimension, so (h, w) is a clean + # multiple of the inner chunk — the Zarr v3 shard-divisibility requirement. + inner_chunks = ( + calculate_aligned_chunk_size(h, 512), + calculate_aligned_chunk_size(w, 512), + ) arr = conditions.create_array( array_name, shape=(h, w), dtype="float32", - chunks=( - calculate_aligned_chunk_size(h, 512), - calculate_aligned_chunk_size(w, 512), - ), + chunks=inner_chunks, + shards=(h, w), compressors=zarr.codecs.BloscCodec(cname="zstd", clevel=5), fill_value=float("nan"), dimension_names=["y", "x"], diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py index 7fea3229..9c86556f 100644 --- a/tests/test_s1_rtc_ingest.py +++ b/tests/test_s1_rtc_ingest.py @@ -2,6 +2,7 @@ from __future__ import annotations +import os from math import ceil from pathlib import Path from unittest.mock import patch @@ -26,6 +27,7 @@ ingest_s1tiling_conditions, parse_s1tiling_filename, ) +from eopf_geozarr.conversion.utils import calculate_aligned_chunk_size # ============================================================================= # Constants @@ -677,6 +679,55 @@ def test_data_integrity_roundtrip( actual = root["ascending"]["conditions"]["gamma_area_037"][:] np.testing.assert_allclose(actual, expected, rtol=1e-6) + def test_gamma_area_is_sharded( + self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path + ) -> None: + """The condition array carries a sharding codec: one shard over the full (y, x) extent, + 512-aligned inner chunks (the same layout vv/vh already use).""" + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=gamma_area_geotiff, + ) + arr = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3)[ + "ascending" + ]["conditions"]["gamma_area_037"] + # shards == full extent (None would mean unsharded — the pre-fix layout) + assert arr.shards == (SIZE, SIZE) + assert arr.chunks == (calculate_aligned_chunk_size(SIZE, 512),) * 2 + + def test_sharding_collapses_chunk_objects_to_one(self, s1_store_with_acquisition: Path) -> None: + """A multi-chunk condition array lands as a SINGLE on-disk shard object, not one object per + inner chunk — the object-count collapse (real gamma_area: ~900 chunk objects → 1 shard).""" + # 1098 sq with a 366 sq inner chunk = 3x3 = 9 inner chunks that, sharded, share one shard. + big = 1098 + rng = np.random.default_rng(7) + data = rng.uniform(0.5, 2.0, (big, big)).astype(np.float32) + gpath = s1_store_with_acquisition.parent / "GAMMA_AREA_BIG_037.tif" + _create_synthetic_geotiff( + gpath, data, transform=from_bounds(XMIN, YMIN, XMAX, YMAX, big, big) + ) + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=gpath, + ) + arr = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3)[ + "ascending" + ]["conditions"]["gamma_area_037"] + assert arr.chunks == (366, 366) + assert arr.shards == (big, big) + # exactly one chunk-data object on disk (the shard), regardless of the 9 inner chunks + array_dir = s1_store_with_acquisition / "ascending" / "conditions" / "gamma_area_037" + data_objects = [ + f for _r, _d, files in os.walk(array_dir) for f in files if f != "zarr.json" + ] + assert len(data_objects) == 1, data_objects + # values still byte-identical through the shard + np.testing.assert_allclose(arr[:], data, rtol=1e-6) + def test_multiple_conditions( self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path, lia_geotiff: Path ) -> None: From 27fe68691e4f80055130e2d4b1260f72d9ab79b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Fri, 19 Jun 2026 19:55:39 +0100 Subject: [PATCH 13/19] fix(s1-rtc): heal a multiscale level missing `time` on append (robust writer) (#200) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * perf(s1-rtc): shard the conditions arrays like the vv/vh pyramid A real S1 RTC cube is 3807 objects / 3.5 GB, of which 3604 (94.7%) are the conditions/gamma_area_ arrays: [10980,10980] float32, inner chunk 366², NO sharding_indexed -> ~900 tiny chunk objects each. They are time-invariant yet dominate the object count, which dominated the ingest's S3 upload wall-time (a live pod sat ~34 min in "Uploading store" at 9 millicores). The vv/vh/border_mask display pyramid is already sharded (one shard per time slice over the full spatial extent, inner 366²). Apply that same existing layout to the condition arrays: add shards=(h, w) to the one create_array in ingest_s1tiling_conditions. All condition arrays (gamma_area, lia, incidence_angle) share that write path and the same 2D full-resolution shape, so all collapse from ~900 chunk objects to 1 shard object (cube ~3807 -> ~210). calculate_aligned_chunk_size returns a divisor of the dimension, so (h, w) is a clean multiple of the inner chunk (Zarr v3 shard-divisibility). conditions arrays are NOT in TiTiler's render path (vv/vh/border_mask), so this does not touch the web-render layout; it only makes a client read a condition array in one ranged GET instead of ~900. Values are byte-identical. Tests: +2 (sharding codec present; 9 inner chunks -> 1 on-disk shard object + byte-identical roundtrip). 57 passed. Spec: claude-docs/specs/s1_gamma_area_sharding.md. Cross-repo Task T5 of data-pipeline/claude-docs/plans/s1_ingest_upload_perf.md. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP * docs(s1-rtc): record real-S3 sharding benchmark in the T5 spec Validated on the live OVH bucket (laptop->DE): object collapse 100->1 (prod ~900->1), PUT 1.7x faster even with batched concurrency on, divisibility valid at the production 10980² (aligned 366 divides 10980), reads byte-identical. Honest caveat recorded: a full-array sequential read is NOT faster sharded (same bytes, one un-parallelizable object) — the win is object-count (upload + listing) and windowed/partial cloud reads, not full-read throughput. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP * refactor(s1-rtc): trim the conditions-sharding comment to match vv/vh style Comment-only. The surrounding vv/vh sharding has no inline explainer; the long cloud-access rationale now lives in claude-docs/specs/s1_gamma_area_sharding.md. Keep only the non-obvious bits: why one shard, and the Zarr v3 shard-divisibility invariant. No behavior change (20 condition/shard tests green). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP * chore(s1-rtc): drop the claude-docs spec from this PR The data-model repo has no claude-docs/specs convention; the spec was noise for this PR's reviewers. The problem statement, real-S3 benchmark and migration note live in the data-pipeline plan + tracking issue EOPF-Explorer/data-pipeline#288 and PR #197's description. Also drop the now-dangling spec path from the code comment (the rationale stays inline). No behavior change (20 condition/shard tests green). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011LsWkVvRfkRzjqAMrzfmRP * fix(s1-rtc): heal a multiscale level missing `time` on append (robust writer) `ingest_s1tiling_acquisition` resized `level["time"]` on every multiscale level, assuming a fresh build created `time` at each level (#192). A cube built before #192 -- or left half-built by an interrupted append -- can carry `r10m/time` yet lack it at a coarser level, so the resize raised `KeyError: 'time'` and, because the append consistency check validated only CRS + shape, the ingest was non-convergent (observed on 30TWM). Before the per-level write loop, recreate any missing-level `time` from `r10m/time` (backfilling the existing slices so prior timestamps are preserved), or raise a clear error when the cube is inconsistent in a way a backfill cannot fix (a level's length disagrees with `r10m/time`, or `r10m` has slices but no `time`). This is the durable upstream counterpart to the data-pipeline guard (data-pipeline #294), making that orchestration-side mitigation belt-and-suspenders. Tests: 4 new cases (heal, no-op when healthy, raise on half-built, raise on missing r10m/time); full s1_ingest suite 61 passed. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01V3qS75byrUuCSHFqcWi26B --------- Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/conversion/s1_ingest.py | 32 ++++++++ tests/test_s1_rtc_ingest.py | 95 ++++++++++++++++++++++++ 2 files changed, 127 insertions(+) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 62d723e0..77b6b0cd 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -684,6 +684,38 @@ def ingest_s1tiling_acquisition( dt_ns = np.datetime64(meta.datetime).astype("datetime64[ns]").astype(np.int64) + # Heal a multiscale level missing `time` before the per-level resize below. A cube built before + # #192 -- or left half-built by an interrupted append -- can carry `r10m/time` yet lack it at a + # coarser level; the unconditional `level["time"].resize` then raises `KeyError: 'time'` and, + # because the consistency check above validates only CRS + shape, the append is non-convergent. + # Recreate the missing-level coordinate from `r10m/time` (backfilling the existing slices so prior + # timestamps are preserved), or raise if the cube is inconsistent in a way a backfill cannot fix. + if "time" in r10m: + ref_time = np.asarray(r10m["time"][:]) + healed = [] + for level_name in data_by_level: + level = orbit[level_name] + if level_name == "r10m" or "time" in level: + continue + level_len = level["vv"].shape[0] + if level_len != ref_time.shape[0]: + raise ValueError( + f"Cannot append to {orbit_direction}/{level_name}: it has {level_len} slice(s) " + f"but r10m/time has {ref_time.shape[0]}; the cube is half-built and `time` cannot " + "be safely backfilled (wipe + reingest)" + ) + _create_time_coordinate_array(level) + level["time"].resize((ref_time.shape[0],)) + level["time"][:] = ref_time + healed.append(level_name) + if healed: + log.info("Healed missing per-level `time`", levels=healed) + elif current_size > 0: + raise ValueError( + f"Cannot append to {orbit_direction}: r10m has {current_size} slice(s) but no `time` " + "coordinate (no backfill source -- wipe + reingest)" + ) + # Write data + the `time` coordinate at all levels (time is replicated per level so datetime # `.sel` resolves at any rendered scale, #192). for level_name, (vv_lev, vh_lev, mask_lev) in data_by_level.items(): diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py index 9c86556f..9471f70b 100644 --- a/tests/test_s1_rtc_ingest.py +++ b/tests/test_s1_rtc_ingest.py @@ -993,3 +993,98 @@ def test_r10m_time_still_int64_for_register( arr = root["ascending"]["r10m"]["time"] assert arr.dtype == np.dtype("int64") assert str(np.datetime64(int(arr[0]), "ns")).startswith("2023-01-15") + + +# ============================================================================= +# Per-level `time` self-heal on append (robust to a pre-#192 / half-built cube) +# ============================================================================= + + +class TestPerLevelTimeHeal: + """The append recreates a multiscale level's missing `time` from r10m/time instead of raising + `KeyError: 'time'` (a cube built before #192, or left half-built by an interrupted append), and + refuses to mis-heal a genuinely inconsistent cube.""" + + def _paths(self, d: Path, stamp: str) -> tuple[Path, Path, Path]: + return ( + d / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC.tif", + d / f"s1a_32TQM_vh_ASC_037_{stamp}_GammaNaughtRTC.tif", + d / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC_BorderMask.tif", + ) + + def _coarse_levels(self, store_path: Path) -> list[str]: + root = zarr.open_group(str(store_path), mode="r", zarr_format=3) + return [n for n, _ in root["ascending"].groups() if n not in ("r10m", "conditions")] + + def _drop_time(self, store_path: Path, level: str) -> None: + """Simulate a pre-#192 cube by removing a level's `time` array. `ingest_s1tiling_acquisition` + does not consolidate, so a filesystem removal is enough for the group to no longer see it.""" + import shutil + + shutil.rmtree(Path(store_path) / "ascending" / level / "time") + + def test_append_heals_levels_missing_time( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """A coarser level lacking `time` is recreated from r10m/time (prior slices preserved); the + append that previously crashed now succeeds.""" + a1 = self._paths(s1_geotiff_dir, "20230115t061234") + a2 = self._paths(s1_geotiff_dir, "20230127t061235") + ingest_s1tiling_acquisition(*a1, s1_store_path, "ascending") + coarse = self._coarse_levels(s1_store_path) + assert coarse # sanity: there ARE coarser levels + for lvl in coarse: + self._drop_time(s1_store_path, lvl) + + idx = ingest_s1tiling_acquisition(*a2, s1_store_path, "ascending") # was KeyError: 'time' + + assert idx == 1 + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + ref = list(root["ascending"]["r10m"]["time"][:]) + assert len(ref) == 2 + for lvl in coarse: + t = root["ascending"][lvl]["time"] + assert t.dtype == np.dtype("int64") + assert list(t[:]) == ref # backfilled prior slice + appended new slice, matching r10m + + def test_append_noop_when_all_levels_have_time( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """A healthy cube: the heal is a no-op and the append works normally.""" + a1 = self._paths(s1_geotiff_dir, "20230115t061234") + a2 = self._paths(s1_geotiff_dir, "20230127t061235") + ingest_s1tiling_acquisition(*a1, s1_store_path, "ascending") + idx = ingest_s1tiling_acquisition(*a2, s1_store_path, "ascending") + assert idx == 1 + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + for lvl in self._coarse_levels(s1_store_path): + assert root["ascending"][lvl]["time"].shape[0] == 2 + + def test_append_raises_on_half_built_cube( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """A level whose data length disagrees with r10m/time (and lacks `time`) is unhealable -> raise + rather than write a wrong-length coordinate.""" + a1 = self._paths(s1_geotiff_dir, "20230115t061234") + a2 = self._paths(s1_geotiff_dir, "20230127t061235") + ingest_s1tiling_acquisition(*a1, s1_store_path, "ascending") + ingest_s1tiling_acquisition(*a2, s1_store_path, "ascending") # 2 slices + + root = zarr.open_group(str(s1_store_path), mode="r+", zarr_format=3) + r20m = root["ascending"]["r20m"] + _, h, w = r20m["vv"].shape + r20m["vv"].resize((1, h, w)) # half-built: r20m has 1 slice, r10m has 2 + self._drop_time(s1_store_path, "r20m") + + with pytest.raises(ValueError, match="half-built"): + ingest_s1tiling_acquisition(*a1, s1_store_path, "ascending") + + def test_append_raises_when_r10m_time_missing( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """r10m holds slices but no `time` -> no backfill source -> raise (not a silent invention).""" + a1 = self._paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(*a1, s1_store_path, "ascending") + self._drop_time(s1_store_path, "r10m") + with pytest.raises(ValueError, match="no backfill source"): + ingest_s1tiling_acquisition(*a1, s1_store_path, "ascending") From c6baebfae568a1f5538636a175048d8ddffcb8e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Mon, 22 Jun 2026 17:21:59 +0100 Subject: [PATCH 14/19] fix(s1-rtc): write CF `_FillValue` on float arrays (#172 parity) (#201) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * fix(s1-rtc): write CF `_FillValue` on float arrays (data-model #172 parity) The S1 RTC ingest writer (`s1_ingest.py`) was never touched by #172, which added the CF `_FillValue` attribute to the S2 and S1-GRD paths so xarray can mask NaN nodata via `use_zarr_fill_value_as_mask=True` despite xarray #11345 — the zarr-level `fill_value` field alone is not surfaced through xarray's encoding. Consequently S1 RTC cubes carried only `grid_mapping` on their float bands, while S2 cubes carry `_FillValue` (`AAAAAAAA+H8=`) + standard_name + units, leaving S1 nodata unmaskable for generic xarray/CF consumers. Set the CF `_FillValue` (FillValueCoder-encoded, matching S2) plus standard_name and units on the float backscatter bands (vv/vh) at every multiscale level, and `_FillValue` on the float condition arrays (gamma_area, lia). The two near-identical vv/vh creation loops — new store vs. new orbit on an existing store — are unified into a shared `_create_band_arrays` helper so they can't drift again (the inline path previously lacked the metadata). Not ported (S2-specific or N/A to the RTC writer): the scale/offset + CastValue codec (S1 vv/vh are float32, unpacked), v0 deprecation (RTC is v1-only), and `sanitize_array_attrs` (s1_ingest writes controlled attrs zarr-direct and never copies source EOPF attrs). Tests: create-store unit, end-to-end ingest across both orbit-creation paths, and float-conditions coverage. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011sBdb17MWSAn7VEPSnqKJU * refactor(s1-rtc): address review — round-trip test, orbit-builder dedup, typed attrs - Add an end-to-end NaN-masking round-trip test (open with `use_zarr_fill_value_as_mask=True`, `to_masked_array()`), mirroring the S2 guarantee in `test_array_attrs.py`. Proves `_FillValue` actually masks nodata, not just that the attribute is present. - Unify the two orbit-creation paths — new store (`create_s1_store`) and new orbit added to an existing store (`ingest_s1tiling_acquisition`) — into a shared `_build_orbit_group` helper. Removes ~60 lines of duplication and fixes the inline path's missing `proj:code` on level groups (guarded by a new test). Net: `s1_ingest.py` shrinks while gaining the fix. - Model the backscatter band CF attrs as a `S1BackscatterAttrsJSON` TypedDict (PR review request), alongside the existing `Standard*CoordAttrsJSON` types. - Document the internal `xarray.backends.zarr.FillValueCoder` import. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_011sBdb17MWSAn7VEPSnqKJU --------- Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/conversion/s1_ingest.py | 192 ++++++++++----------- src/eopf_geozarr/data_api/geozarr/types.py | 12 ++ tests/test_s1_rtc_ingest.py | 129 ++++++++++++++ 3 files changed, 229 insertions(+), 104 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 77b6b0cd..f8175fbb 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -19,6 +19,7 @@ from dataclasses import dataclass from math import ceil from pathlib import Path +from typing import TYPE_CHECKING, cast import numpy as np import rasterio @@ -26,12 +27,20 @@ import zarr import zarr.codecs from pyproj import CRS as PyprojCRS + +# `FillValueCoder` is xarray's internal CF fill-value encoder: the canonical way to produce +# the base64 `_FillValue` xarray reads back under `use_zarr_fill_value_as_mask` (xarray +# #11345); the S2 path relies on the same mechanism. Internal API — revisit if xarray moves it. +from xarray.backends.zarr import FillValueCoder from zarr_cm import geo_proj from zarr_cm import multiscales as multiscales_cm from zarr_cm import spatial as spatial_cm from eopf_geozarr.conversion.utils import calculate_aligned_chunk_size +if TYPE_CHECKING: + from eopf_geozarr.data_api.geozarr.types import S1BackscatterAttrsJSON + log = structlog.get_logger() # ============================================================================= @@ -308,6 +317,21 @@ def _create_spatial_coordinate_arrays( "standard_name": "time", } +# CF `_FillValue` for the float32 arrays, mirroring the S1 GRD converter (geozarr.py) +# and S2 (data-model #172). This is what lets xarray mask NaN nodata via +# `use_zarr_fill_value_as_mask=True` despite xarray #11345 — the zarr-level `fill_value` +# field alone is not surfaced through xarray's encoding. We encode it with +# `FillValueCoder` so the stored attribute matches the base64 form S2 writes +# (`AAAAAAAA+H8=`). The S2 path gets this for free via `to_zarr`; this store is written +# zarr-direct, so the attribute must be set explicitly. +FLOAT32_NAN_FILL_VALUE = FillValueCoder.encode(np.nan, np.dtype("float32")) +# CF metadata for the backscatter bands (vv/vh). +BACKSCATTER_CF_ATTRS: S1BackscatterAttrsJSON = { + "standard_name": "surface_backwards_scattering_coefficient_of_radar_wave", + "units": "1", + "_FillValue": FLOAT32_NAN_FILL_VALUE, +} + def _create_time_coordinate_array(level_group: zarr.Group) -> None: """Create the CF-encoded ``time`` coordinate (length 0, grown on append) on one level group. @@ -358,20 +382,55 @@ def _add_grid_mapping(group: zarr.Group, crs_string: str) -> None: arr.attrs["grid_mapping"] = "spatial_ref" -def create_s1_store( - store_path: str | Path, - orbit_direction: str, - metadata: S1TilingMetadata, +def _create_band_arrays(level_group: zarr.Group, level_h: int, level_w: int) -> None: + """Create the (time, y, x) data bands (vv, vh, border_mask) for one multiscale level. + + Float backscatter bands (vv/vh) get the CF ``_FillValue``/``standard_name``/``units`` + attributes so readers can mask NaN nodata (xarray #11345) — parity with the S1 GRD + converter and S2 (data-model #172). Shared by ``create_s1_store`` (new store) and + ``ingest_s1tiling_acquisition`` (new orbit added to an existing store) so the two + creation paths cannot drift. + """ + inner_chunks = ( + 1, + calculate_aligned_chunk_size(level_h, 512), + calculate_aligned_chunk_size(level_w, 512), + ) + shard_shape = (1, level_h, level_w) + for name, dtype, fill in [ + ("vv", "float32", float("nan")), + ("vh", "float32", float("nan")), + ("border_mask", "uint8", 0), + ]: + band = level_group.create_array( + name, + shape=(0, level_h, level_w), + dtype=dtype, + chunks=inner_chunks, + shards=shard_shape, + compressors=zarr.codecs.BloscCodec(cname="zstd", clevel=5), + fill_value=fill, + dimension_names=["time", "y", "x"], + ) + if name in ("vv", "vh"): + # cast: zarr's `attrs.update` is typed for its JSON union, which a TypedDict + # (invariant) doesn't satisfy structurally; the values are JSON-safe. + band.attrs.update(cast("dict", BACKSCATTER_CF_ATTRS)) + + +def _build_orbit_group( + root: zarr.Group, orbit_direction: str, metadata: S1TilingMetadata ) -> zarr.Group: - """Create a new S1 GRD RTC Zarr V3 store with full conventions metadata. + """Create one orbit group (asc/desc) with full GeoZarr metadata. - Returns the root group. + Builds the multiscale level groups (bands + spatial coordinates + grid-mapping + + per-level ``time``) and the native-resolution per-acquisition metadata coordinates. + Shared by ``create_s1_store`` (new store) and ``ingest_s1tiling_acquisition`` (new + orbit added to an existing store) so the two creation paths cannot drift — the inline + path previously omitted ``proj:code`` on level groups. """ layout = compute_multiscales_layout(metadata.shape, metadata.spatial_transform) - - root = zarr.open_group(str(store_path), mode="w-", zarr_format=3) orbit_group = root.create_group(orbit_direction) - orbit_group.attrs.update( { "zarr_conventions": ZARR_CONVENTIONS, @@ -398,28 +457,7 @@ def create_s1_store( } ) - inner_chunks = ( - 1, - calculate_aligned_chunk_size(level_h, 512), - calculate_aligned_chunk_size(level_w, 512), - ) - shard_shape = (1, level_h, level_w) - - for name, dtype, fill in [ - ("vv", "float32", float("nan")), - ("vh", "float32", float("nan")), - ("border_mask", "uint8", 0), - ]: - level_group.create_array( - name, - shape=(0, level_h, level_w), - dtype=dtype, - chunks=inner_chunks, - shards=shard_shape, - compressors=zarr.codecs.BloscCodec(cname="zstd", clevel=5), - fill_value=fill, - dimension_names=["time", "y", "x"], - ) + _create_band_arrays(level_group, level_h, level_w) _create_spatial_coordinate_arrays( level_group, level_h, level_w, level_entry["spatial:transform"] @@ -450,6 +488,20 @@ def create_s1_store( fill_value="", dimension_names=["time"], ) + return orbit_group + + +def create_s1_store( + store_path: str | Path, + orbit_direction: str, + metadata: S1TilingMetadata, +) -> zarr.Group: + """Create a new S1 GRD RTC Zarr V3 store with full conventions metadata. + + Returns the root group. + """ + root = zarr.open_group(str(store_path), mode="w-", zarr_format=3) + _build_orbit_group(root, orbit_direction, metadata) log.info( "Created S1 store", @@ -559,79 +611,9 @@ def ingest_s1tiling_acquisition( else: root = zarr.open_group(str(store_path), mode="r+", zarr_format=3) if orbit_direction not in root: - # Create orbit direction group in existing store - orbit_group = root.create_group(orbit_direction) - layout = compute_multiscales_layout(meta.shape, meta.spatial_transform) - orbit_group.attrs.update( - { - "zarr_conventions": ZARR_CONVENTIONS, - "multiscales": { - "layout": layout, - "resampling_method": "average", - }, - "proj:code": meta.crs, - "spatial:dimensions": ["y", "x"], - "spatial:bbox": meta.bounds, - } - ) - for level_entry in layout: - level_name = level_entry["asset"] - level_h, level_w = level_entry["spatial:shape"] - level_group = orbit_group.create_group(level_name) - level_group.attrs.update( - { - "spatial:shape": [level_h, level_w], - "spatial:transform": level_entry["spatial:transform"], - } - ) - inner_chunks = ( - 1, - calculate_aligned_chunk_size(level_h, 512), - calculate_aligned_chunk_size(level_w, 512), - ) - shard_shape = (1, level_h, level_w) - for name, dtype, fill in [ - ("vv", "float32", float("nan")), - ("vh", "float32", float("nan")), - ("border_mask", "uint8", 0), - ]: - level_group.create_array( - name, - shape=(0, level_h, level_w), - dtype=dtype, - chunks=inner_chunks, - shards=shard_shape, - compressors=zarr.codecs.BloscCodec(cname="zstd", clevel=5), - fill_value=fill, - dimension_names=["time", "y", "x"], - ) - _create_spatial_coordinate_arrays( - level_group, level_h, level_w, level_entry["spatial:transform"] - ) - _add_grid_mapping(level_group, meta.crs) - # `time` coordinate on every level (datetime `.sel` at any scale, #192). - _create_time_coordinate_array(level_group) - r10m = orbit_group["r10m"] - for name, dtype, fill in [ - ("absolute_orbit", "int32", 0), - ("relative_orbit", "int32", 0), - ]: - r10m.create_array( - name, - shape=(0,), - dtype=dtype, - chunks=(512,), - fill_value=fill, - dimension_names=["time"], - ) - r10m.create_array( - "platform", - shape=(0,), - dtype=" None: + """Float backscatter bands must declare a CF ``_FillValue`` attribute at + *every* multiscale level, matching S2 (data-model #172 / xarray #11345). + + The zarr-level ``fill_value`` alone is not surfaced by xarray's encoding, + so ``to_masked_array()`` / ``use_zarr_fill_value_as_mask=True`` cannot mask + NaN nodata without the attribute. ``test_array_attrs`` only guards the S2 + ``/measurements/`` layout, so the S1 RTC layout was previously unchecked. + """ + from xarray.backends.zarr import FillValueCoder + + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + orbit = root["ascending"] + expected = FillValueCoder.encode(np.nan, np.dtype("float32")) + for level_name, _, _ in OVERVIEW_CHAIN: + for band in ("vv", "vh"): + attrs = dict(orbit[level_name][band].attrs) + assert attrs.get("_FillValue") == expected, ( + f"{level_name}/{band} missing/!= CF _FillValue" + ) + assert ( + attrs.get("standard_name") + == "surface_backwards_scattering_coefficient_of_radar_wave" + ) + assert attrs.get("units") == "1" + def test_no_tile_matrix_set( self, s1_store_path: Path, sample_metadata: S1TilingMetadata ) -> None: @@ -354,6 +382,83 @@ def test_second_acquisition_appends(self, s1_geotiff_dir: Path, s1_store_path: P root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) assert root["ascending"]["r10m"]["vv"].shape[0] == 2 + def test_ingested_bands_declare_cf_fill_value( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """End-to-end (production path): vv/vh carry CF ``_FillValue``/``standard_name``/ + ``units`` at every level for BOTH the store-creating orbit (``create_s1_store``) + and a second orbit added via the inline new-orbit path — the two paths that were + previously inconsistent. Parity with S2 / S1 GRD (#172; xarray #11345). + """ + from xarray.backends.zarr import FillValueCoder + + vv, vh, mask = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") # create_s1_store + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "descending") # inline new orbit + expected = FillValueCoder.encode(np.nan, np.dtype("float32")) + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + for orbit in ("ascending", "descending"): + for level_name, _, _ in OVERVIEW_CHAIN: + for band in ("vv", "vh"): + attrs = dict(root[orbit][level_name][band].attrs) + assert attrs.get("_FillValue") == expected, f"{orbit}/{level_name}/{band}" + assert ( + attrs.get("standard_name") + == "surface_backwards_scattering_coefficient_of_radar_wave" + ) + assert attrs.get("units") == "1" + + def test_new_orbit_level_groups_carry_proj_code( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """A second orbit added to an existing store must get the same per-level metadata + as the store-creating orbit — incl. ``proj:code`` on every level group. The inline + new-orbit path previously omitted it (drift vs ``create_s1_store``).""" + vv, vh, mask = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "descending") + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + for orbit in ("ascending", "descending"): + for level_name, _, _ in OVERVIEW_CHAIN: + attrs = dict(root[orbit][level_name].attrs) + assert attrs.get("proj:code") == CRS, f"{orbit}/{level_name} missing proj:code" + + def test_fill_value_masking_roundtrip(self, tmp_path: Path, s1_store_path: Path) -> None: + """End-to-end: NaN nodata in vv comes back masked when the cube is reopened with + ``use_zarr_fill_value_as_mask=True`` — the behaviour the CF ``_FillValue`` attribute + exists to enable despite xarray #11345. Mirrors the S2 guarantee in + ``tests/test_array_attrs.py::test_fill_value_masking_roundtrip``. + """ + stamp = "20230115t061234" + rng = np.random.default_rng(0) + vv_data = rng.uniform(0.0, 1.0, (SIZE, SIZE)).astype(np.float32) + vv_data[0:16, 0:16] = np.nan # nodata patch + vh_data = rng.uniform(0.0, 0.5, (SIZE, SIZE)).astype(np.float32) + mask_data = np.ones((SIZE, SIZE), dtype=np.uint8) + vv = tmp_path / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC.tif" + vh = tmp_path / f"s1a_32TQM_vh_ASC_037_{stamp}_GammaNaughtRTC.tif" + mask = tmp_path / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC_BorderMask.tif" + _create_synthetic_geotiff(vv, vv_data, tags=ACQ1_TAGS) + _create_synthetic_geotiff(vh, vh_data, tags=ACQ1_TAGS) + _create_synthetic_geotiff(mask, mask_data, tags=ACQ1_TAGS) + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") + + ds = xr.open_dataset( + str(s1_store_path / "ascending" / "r10m"), + engine="zarr", + consolidated=False, + decode_times=False, + decode_coords=False, + use_zarr_fill_value_as_mask=True, + ) + try: + masked = ds["vv"].to_masked_array() + assert np.ma.is_masked(masked), "NaN nodata must be masked via `_FillValue`" + assert masked.mask[0, 0, 0], "nodata cell must be masked" + assert not masked.mask[0, -1, -1], "valid cell must not be masked" + finally: + ds.close() + def test_preserves_data_integrity(self, s1_geotiff_dir: Path, s1_store_path: Path) -> None: vv, vh, mask = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") @@ -679,6 +784,30 @@ def test_data_integrity_roundtrip( actual = root["ascending"]["conditions"]["gamma_area_037"][:] np.testing.assert_allclose(actual, expected, rtol=1e-6) + def test_float_conditions_declare_cf_fill_value( + self, + s1_store_with_acquisition: Path, + gamma_area_geotiff: Path, + lia_geotiff: Path, + ) -> None: + """Float condition arrays (gamma_area, lia) must declare a CF ``_FillValue`` so + readers mask NaN nodata (xarray #11345), like vv/vh (#172).""" + from xarray.backends.zarr import FillValueCoder + + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=gamma_area_geotiff, + lia_path=lia_geotiff, + ) + expected = FillValueCoder.encode(np.nan, np.dtype("float32")) + conditions = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3)[ + "ascending" + ]["conditions"] + for arr_name in ("gamma_area_037", "lia_037"): + assert dict(conditions[arr_name].attrs).get("_FillValue") == expected, arr_name + def test_gamma_area_is_sharded( self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path ) -> None: From 9bd5a551a63bd5872823cdd98c4db60958995f74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Tue, 23 Jun 2026 11:42:05 +0100 Subject: [PATCH 15/19] fix(s1-rtc): store nodata as NaN, not 0, so titiler masks it transparent (#202) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit S1 RTC previews render the out-of-swath region as opaque black while the S2 reference renders it transparent. The cause is the *data*, not the metadata: #201 already gave vv/vh the same dtype/`fill_value`=NaN/ `_FillValue`/`grid_mapping` encoding as S2, but those are inert because no NaN is ever written — s1tiling stores `0.0` out of swath, and titiler treats `0` as valid data and paints it black. Mask nodata to NaN at the writer: - vv/vh: `np.where(border_mask == 0, NaN, ...)` — border_mask is the authoritative valid-data mask (0 = no-data). `_downsample_2d` already uses `np.nanmean` for floats, so NaN propagates to every overview level. - float conditions (gamma_area/lia): masked read off the GeoTIFF's declared nodata (border_mask is N/A for static conditions); a no-op when no nodata is declared. Only newly re-ingested cubes get NaN; existing cubes are remediated separately. Tests: NaN ⟺ border_mask==0 at native + overview levels, masking round-trips via `use_zarr_fill_value_as_mask`, conditions declared-nodata → NaN. Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/conversion/s1_ingest.py | 12 ++- tests/test_s1_rtc_ingest.py | 99 ++++++++++++++++++++---- 2 files changed, 97 insertions(+), 14 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index f8175fbb..5feed9e2 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -646,6 +646,13 @@ def ingest_s1tiling_acquisition( vv_max=float(np.nanmax(vv_data)), ) + # nodata → NaN: s1tiling writes 0 out of swath, which titiler treats as valid data and renders + # opaque black. Store NaN there instead so it masks transparent like the S2 reference. The + # border mask is the authoritative valid-data mask (0 = no-data); `_downsample_2d` uses + # `np.nanmean` for floats, so NaN propagates to every overview level below. + vv_data = np.where(mask_data == 0, np.nan, vv_data).astype("float32") + vh_data = np.where(mask_data == 0, np.nan, vh_data).astype("float32") + # Determine time index r10m = orbit["r10m"] current_size = r10m["vv"].shape[0] @@ -990,7 +997,10 @@ def ingest_s1tiling_conditions( array_name = f"{label}_{orbit_suffix}" with _rasterio_env(cond_path), rasterio.open(str(cond_path)) as src: - data = src.read(1).astype(np.float32) + # nodata → NaN via the GeoTIFF's declared nodata (border_mask is N/A for static + # conditions), so out-of-coverage pixels mask transparent like vv/vh. A no-op when + # the GeoTIFF declares no nodata. + data = src.read(1, masked=True).filled(np.nan).astype(np.float32) h, w = data.shape diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py index 6c7cc1dc..c70addf0 100644 --- a/tests/test_s1_rtc_ingest.py +++ b/tests/test_s1_rtc_ingest.py @@ -68,8 +68,9 @@ def _create_synthetic_geotiff( crs: str = CRS, transform: rasterio.transform.Affine | None = None, tags: dict[str, str] | None = None, + nodata: float | None = None, ) -> None: - """Write a single-band GeoTIFF with optional metadata tags.""" + """Write a single-band GeoTIFF with optional metadata tags and declared nodata.""" if transform is None: transform = TRANSFORM with rasterio.open( @@ -82,6 +83,7 @@ def _create_synthetic_geotiff( dtype=data.dtype, crs=crs, transform=transform, + nodata=nodata, ) as dst: if tags: dst.update_tags(**tags) @@ -424,17 +426,22 @@ def test_new_orbit_level_groups_carry_proj_code( assert attrs.get("proj:code") == CRS, f"{orbit}/{level_name} missing proj:code" def test_fill_value_masking_roundtrip(self, tmp_path: Path, s1_store_path: Path) -> None: - """End-to-end: NaN nodata in vv comes back masked when the cube is reopened with - ``use_zarr_fill_value_as_mask=True`` — the behaviour the CF ``_FillValue`` attribute - exists to enable despite xarray #11345. Mirrors the S2 guarantee in + """End-to-end: out-of-swath nodata (``border_mask == 0``) comes back masked when the cube + is reopened with ``use_zarr_fill_value_as_mask=True`` — the behaviour the CF ``_FillValue`` + attribute exists to enable despite xarray #11345. Mirrors the S2 guarantee in ``tests/test_array_attrs.py::test_fill_value_masking_roundtrip``. + + The nodata region comes from ``border_mask``, not a pre-seeded NaN: s1tiling stores ``0.0`` + out of swath, and the writer is what must convert that to NaN. """ stamp = "20230115t061234" rng = np.random.default_rng(0) - vv_data = rng.uniform(0.0, 1.0, (SIZE, SIZE)).astype(np.float32) - vv_data[0:16, 0:16] = np.nan # nodata patch - vh_data = rng.uniform(0.0, 0.5, (SIZE, SIZE)).astype(np.float32) + vv_data = rng.uniform(0.1, 1.0, (SIZE, SIZE)).astype(np.float32) + vh_data = rng.uniform(0.1, 0.5, (SIZE, SIZE)).astype(np.float32) mask_data = np.ones((SIZE, SIZE), dtype=np.uint8) + mask_data[0:16, 0:16] = 0 # out-of-swath border + vv_data[0:16, 0:16] = 0.0 # s1tiling stores 0 where there is no swath + vh_data[0:16, 0:16] = 0.0 vv = tmp_path / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC.tif" vh = tmp_path / f"s1a_32TQM_vh_ASC_037_{stamp}_GammaNaughtRTC.tif" mask = tmp_path / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC_BorderMask.tif" @@ -453,12 +460,49 @@ def test_fill_value_masking_roundtrip(self, tmp_path: Path, s1_store_path: Path) ) try: masked = ds["vv"].to_masked_array() - assert np.ma.is_masked(masked), "NaN nodata must be masked via `_FillValue`" - assert masked.mask[0, 0, 0], "nodata cell must be masked" + assert np.ma.is_masked(masked), "out-of-swath nodata must be masked via `_FillValue`" + assert masked.mask[0, 0, 0], "nodata cell (border_mask==0) must be masked" assert not masked.mask[0, -1, -1], "valid cell must not be masked" finally: ds.close() + def test_nodata_masked_to_nan(self, tmp_path: Path, s1_store_path: Path) -> None: + """The writer stores NaN — not 0 — wherever ``border_mask == 0``, at the native level and + every overview, so titiler masks out-of-swath nodata transparent like the S2 reference. + + NaN must coincide exactly with ``border_mask == 0``: valid pixels stay finite. Root cause + of the "black area" render bug: s1tiling writes 0 out of swath, and 0 is valid data to + titiler. ``np.nanmean`` downsampling must carry the NaN to every overview level. + """ + stamp = "20230115t061234" + rng = np.random.default_rng(7) + vv_data = rng.uniform(0.1, 1.0, (SIZE, SIZE)).astype(np.float32) + vh_data = rng.uniform(0.1, 0.5, (SIZE, SIZE)).astype(np.float32) + mask_data = np.ones((SIZE, SIZE), dtype=np.uint8) + mask_data[0:32, :] = 0 # out-of-swath border band (whole rows) + vv_data[0:32, :] = 0.0 # s1tiling stores 0 there — valid data to titiler today + vh_data[0:32, :] = 0.0 + vv = tmp_path / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC.tif" + vh = tmp_path / f"s1a_32TQM_vh_ASC_037_{stamp}_GammaNaughtRTC.tif" + mask = tmp_path / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC_BorderMask.tif" + _create_synthetic_geotiff(vv, vv_data, tags=ACQ1_TAGS) + _create_synthetic_geotiff(vh, vh_data, tags=ACQ1_TAGS) + _create_synthetic_geotiff(mask, mask_data, tags=ACQ1_TAGS) + ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") + + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + nodata = mask_data == 0 + for band in ("vv", "vh"): + native = root["ascending"]["r10m"][band][0, :, :] + assert np.all(np.isnan(native[nodata])), f"{band}: nodata region must be NaN, not 0" + assert not np.any(np.isnan(native[~nodata])), f"{band}: valid region must stay finite" + assert not np.any(native[nodata] == 0.0), f"{band}: nodata must not read back as 0" + + # NaN propagates through np.nanmean downsampling: the all-nodata top band stays NaN. + coarse = root["ascending"]["r20m"]["vv"][0, :, :] + assert np.isnan(coarse[0, 0]), "all-nodata block must stay NaN at the overview level" + assert not np.any(np.isnan(coarse[-1, :])), "fully-valid bottom row must stay finite" + def test_preserves_data_integrity(self, s1_geotiff_dir: Path, s1_store_path: Path) -> None: vv, vh, mask = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") @@ -466,13 +510,18 @@ def test_preserves_data_integrity(self, s1_geotiff_dir: Path, s1_store_path: Pat # Read back and compare with rasterio.open(str(vv)) as src: expected_vv = src.read(1) + with rasterio.open(str(mask)) as src: + expected_mask = src.read(1).astype(np.uint8) root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) actual_vv = root["ascending"]["r10m"]["vv"][0, :, :] - np.testing.assert_allclose(actual_vv, expected_vv, rtol=1e-6) - # Mask should be exact - with rasterio.open(str(mask)) as src: - expected_mask = src.read(1).astype(np.uint8) + # Valid pixels (border_mask == 1) are preserved exactly; out-of-swath pixels + # (border_mask == 0) are written as NaN, not the raw 0 — the render-bug fix. + valid = expected_mask == 1 + np.testing.assert_allclose(actual_vv[valid], expected_vv[valid], rtol=1e-6) + assert np.all(np.isnan(actual_vv[~valid])), "out-of-swath nodata must be NaN" + + # border_mask itself is stored verbatim (uint8, never masked). actual_mask = root["ascending"]["r10m"]["border_mask"][0, :, :] np.testing.assert_array_equal(actual_mask, expected_mask) @@ -784,6 +833,30 @@ def test_data_integrity_roundtrip( actual = root["ascending"]["conditions"]["gamma_area_037"][:] np.testing.assert_allclose(actual, expected, rtol=1e-6) + def test_conditions_nodata_masked_to_nan( + self, s1_store_with_acquisition: Path, tmp_path: Path + ) -> None: + """A condition GeoTIFF's declared-nodata pixels read back as NaN (not the raw sentinel), + so the auxiliary arrays mask transparent like vv/vh. border_mask is N/A for static + conditions, so the writer relies on the GeoTIFF's declared nodata via a masked read. + """ + data = np.full((SIZE, SIZE), 1.5, dtype=np.float32) + data[0:20, 0:20] = 0.0 # declared-nodata region + cond_path = tmp_path / "GAMMA_AREA_32TQM_037.tif" + _create_synthetic_geotiff(cond_path, data, nodata=0.0) + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=cond_path, + ) + conditions = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3)[ + "ascending" + ]["conditions"] + arr = conditions["gamma_area_037"][:] + assert np.all(np.isnan(arr[0:20, 0:20])), "declared-nodata region must be NaN" + assert not np.any(np.isnan(arr[20:, 20:])), "valid region must stay finite" + def test_float_conditions_declare_cf_fill_value( self, s1_store_with_acquisition: Path, From 79ce15e78a5f07c88cb42f8a9ee08d7ee77302a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Tue, 23 Jun 2026 15:21:51 +0100 Subject: [PATCH 16/19] fix(s1-rtc): consolidate every orbit group, not just the last-ingested one (#203) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit S1 RTC cubes ended up with only one orbit consolidated on disk (staging: 30TWM asc✓/desc✗, 32TNN/30UVU asc✗/desc✓). Root cause: the pipeline ingests acquisitions one orbit at a time in separate pods, each of which strips all consolidated metadata (so `time` can resize), ingests its orbit, then calls `consolidate_s1_store(store, orbit_direction)` — which only re-consolidated the single orbit passed. So whichever orbit was ingested last is the only one left with on-disk consolidated metadata. Consolidate every orbit group present (iterate `root.groups()`), then the root. The minimal append-fetch already pulls all `zarr.json`, so both orbits' metadata is local — same reason the root consolidation already works. Signature unchanged (`orbit_direction` kept for logging/callers). Impact is cosmetic: the root consolidation is complete for both orbits and readers opening at the root get a synthesized child view (titiler renders the unconsolidated orbit fine). The fix matters for clients opening a single orbit group standalone (the STAC `.zarr/` hrefs) with `consolidated=True`, which otherwise fall back to a listing. Test asserts each orbit group is consolidated *standalone* — checking via `root[orbit]` is a false-green because a consolidated root synthesizes the child's view. Existing cubes self-heal on their next re-ingest. Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/conversion/s1_ingest.py | 13 +++++++++++-- tests/test_s1_rtc_ingest.py | 20 ++++++++++++++++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py index 5feed9e2..2b541cae 100644 --- a/src/eopf_geozarr/conversion/s1_ingest.py +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -743,12 +743,21 @@ def ingest_s1tiling_acquisition( def consolidate_s1_store(store_path: str | Path, orbit_direction: str) -> None: - """Consolidate metadata at orbit direction and root levels. + """Consolidate metadata for every orbit-direction group and the root. Must be called AFTER all ingestions complete — consolidated metadata caches array shapes and will become stale if called mid-ingestion. + + Consolidates *every* orbit group present, not just ``orbit_direction``: the + pipeline ingests acquisitions one orbit at a time after stripping all + consolidated metadata (so ``time`` can resize), so consolidating only the + passed orbit would leave the other orbit's group unconsolidated on disk + (readers opening that orbit standalone then fall back to a listing). + ``orbit_direction`` is retained for logging / caller compatibility. """ - zarr.consolidate_metadata(str(store_path), path=orbit_direction, zarr_format=3) + root = zarr.open_group(str(store_path), mode="r", zarr_format=3) + for orbit_name, _ in root.groups(): + zarr.consolidate_metadata(str(store_path), path=orbit_name, zarr_format=3) zarr.consolidate_metadata(str(store_path), zarr_format=3) log.info( "Metadata consolidated", diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py index c70addf0..131939d3 100644 --- a/tests/test_s1_rtc_ingest.py +++ b/tests/test_s1_rtc_ingest.py @@ -652,6 +652,26 @@ def test_consolidate_after_all_ingestions( r10m = root["ascending"]["r10m"] assert r10m["vv"].shape[0] == 2 + def test_consolidate_all_orbits_present( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + """``consolidate_s1_store`` must leave EVERY orbit group consolidated on disk, not just the + one passed. The pipeline ingests acquisitions one orbit at a time after stripping all + consolidated metadata (so ``time`` can resize), so consolidating only the passed orbit left + staging cubes asc✓/desc✗. Each orbit group is checked **standalone**: a consolidated root + synthesises the child's view, so ``root[orbit].metadata.consolidated_metadata`` is a + false-green (non-None even when ``/zarr.json`` lacks it). + """ + vv1, vh1, mask1 = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") + vv2, vh2, mask2 = self._get_acq_paths(s1_geotiff_dir, "20230127t061235") + ingest_s1tiling_acquisition(vv1, vh1, mask1, s1_store_path, "ascending") + ingest_s1tiling_acquisition(vv2, vh2, mask2, s1_store_path, "descending") + consolidate_s1_store(s1_store_path, "descending") # only one orbit passed + + for orbit in ("ascending", "descending"): + grp = zarr.open_group(str(s1_store_path / orbit), mode="r", zarr_format=3) + assert grp.metadata.consolidated_metadata is not None, f"{orbit} orbit not consolidated" + def _get_acq_paths(self, geotiff_dir: Path, stamp: str) -> tuple[Path, Path, Path]: vv = geotiff_dir / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC.tif" vh = geotiff_dir / f"s1a_32TQM_vh_ASC_037_{stamp}_GammaNaughtRTC.tif" From 0734ff1608c3867d4d0a4069926d43c86f0a114a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:42:13 +0100 Subject: [PATCH 17/19] fix(s1-rtc): per-band rescale for the RGB composite render (#204) The render expression emits 3 bands (vv; vh; vv/vh) but supplied a single rescale pair [0,0.2] applied to all three. The vv/vh ratio (natural range ~1-15) saturated under that stretch -> flat blue/purple wash, and low cross-pol water dropped to transparent, which made (correctly geolocated) ascending swaths look mislocated. Give each band its own pair: vv [0,0.4], vh [0,0.1], ratio [1,15]. Cosmetic preview-only change; pixel placement is unchanged. Verified live against the raster API and the S1 STAC tests. Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/stac/s1_rtc.py | 12 ++++++++---- tests/test_s1_stac.py | 2 +- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index 6e50d33e..9e707d87 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -38,16 +38,20 @@ def _rgb_render(orbit: str) -> dict[str, object]: Produces a 3-band false-colour composite (R=VV, G=VH, B=VV/VH ratio) that titiler renders into previews/tiles. ``bidx=[1]`` selects the single time - slice from each multi-band variable; the single ``rescale`` pair (linear - gamma0 units) is applied to every band. + slice from each multi-band variable. Each band gets its own ``rescale`` pair: + VV and VH are linear gamma0 (low values) while the VV/VH ratio spans ~1-15, so + a single shared pair blew the ratio band out to a flat blue/purple wash (and + dropped low-cross-pol water to transparent, making swaths look mislocated). + Per-band stretches keep the composite natural and the swath readable. """ vv = f"/{orbit}:vv" vh = f"/{orbit}:vh" return { "title": "VV, VH, VV/VH composite", "expression": f"{vv};{vh};({vv})/({vh})", - # Linear gamma0 stretch tuned for the S1 RTC product (the 0.0,0.1 default was too bright). - "rescale": [[0.0, 0.2]], + # Per-band linear stretch: VV/VH are low-valued gamma0; the VV/VH ratio spans ~1-15. + # One shared pair saturated the ratio band (purple wash) and dropped low-cross-pol water. + "rescale": [[0.0, 0.4], [0.0, 0.1], [1.0, 15.0]], "bidx": [1], "tilesize": 256, } diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index 732d880f..061f13f8 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -324,7 +324,7 @@ def test_render_extension_rgb_composite(tmp_path: Path) -> None: rgb = item.properties["renders"]["rgb"] assert rgb["expression"] == "/descending:vv;/descending:vh;(/descending:vv)/(/descending:vh)" - assert rgb["rescale"] == [[0.0, 0.2]] + assert rgb["rescale"] == [[0.0, 0.4], [0.0, 0.1], [1.0, 15.0]] assert rgb["bidx"] == [1] assert rgb["tilesize"] == 256 From 2e702ca4f82798a0adb73a250255ba814300d9a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Tue, 23 Jun 2026 20:41:55 +0100 Subject: [PATCH 18/19] feat(s1-rtc): add grid:code (MGRS tile) + grid extension to STAC items (#205) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Both the cube item and (by inheritance) its per-acquisition items now carry the STAC grid extension with grid:code = "MGRS-{tile}". This makes the tile a first-class queryable property: the acquisitions collection becomes natively tile-filterable in STAC Browser, and cube↔acquisition cross-links can filter on grid:code instead of an id-prefix LIKE. Cube-builder only — per-acquisition items inherit it via the existing build_s1_rtc_per_acquisition_items copy ({**base_dict} carries stac_extensions; the per-acq property denylist excludes grid:code). Claude-Session: https://claude.ai/code/session_019z3eVtkSNHhN9vHd8QqGcf Co-authored-by: Claude Opus 4.8 --- src/eopf_geozarr/stac/s1_rtc.py | 6 +++++- tests/test_s1_stac.py | 8 ++++++++ tests/test_s1_stac_per_acquisition.py | 8 ++++++++ 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/eopf_geozarr/stac/s1_rtc.py b/src/eopf_geozarr/stac/s1_rtc.py index 9e707d87..4848eb3e 100644 --- a/src/eopf_geozarr/stac/s1_rtc.py +++ b/src/eopf_geozarr/stac/s1_rtc.py @@ -17,6 +17,7 @@ RENDER_EXT = "https://stac-extensions.github.io/render/v1.0.0/schema.json" DATACUBE_EXT = "https://stac-extensions.github.io/datacube/v2.2.0/schema.json" TIMESTAMPS_EXT = "https://stac-extensions.github.io/timestamps/v1.1.0/schema.json" +GRID_EXT = "https://stac-extensions.github.io/grid/v1.1.0/schema.json" ZARR_MEDIA_TYPE = "application/vnd.zarr; version=3" @@ -216,6 +217,9 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: # Projection extension "proj:code": preferred_proj_code, "proj:bbox": preferred_bbox, + # Grid extension: the Sentinel-2 MGRS tile this cube is gridded onto — a queryable tile id + # (enables tile-filtering the acquisitions collection and cube↔acquisition cross-links). + "grid:code": f"MGRS-{tile_id}", # Render extension: dual-pol RGB composite for previews/tiles (defaults to the preferred orbit) "renders": {"rgb": _rgb_render(preferred_orbit)}, } @@ -224,7 +228,7 @@ def build_s1_rtc_stac_item(zarr_store: str, collection_id: str) -> pystac.Item: if preferred["transform"] is not None: properties["proj:transform"] = preferred["transform"] - stac_extensions = [SAR_EXT, PROJ_EXT, RENDER_EXT, DATACUBE_EXT, TIMESTAMPS_EXT] + stac_extensions = [SAR_EXT, PROJ_EXT, RENDER_EXT, DATACUBE_EXT, TIMESTAMPS_EXT, GRID_EXT] # sat:orbit_state is single-valued, so it's only meaningful when the cube holds a single orbit. A # dual-orbit cube would mislabel half its slices — omit it there (per-acquisition items, which are diff --git a/tests/test_s1_stac.py b/tests/test_s1_stac.py index 061f13f8..99abf6dd 100644 --- a/tests/test_s1_stac.py +++ b/tests/test_s1_stac.py @@ -92,6 +92,14 @@ def test_item_id_matches_tile_id(tmp_path: Path) -> None: assert item.id == "s1-rtc-31TCH" +def test_grid_code_and_extension(tmp_path: Path) -> None: + """The cube declares the grid extension + grid:code = MGRS-{tile} (a queryable tile id).""" + store = _make_s1_store(tmp_path, {"descending": [(T1_NS, "S1A")]}) + item = build_s1_rtc_stac_item(str(store), "sentinel-1-grd-rtc-staging") + assert item.properties["grid:code"] == "MGRS-31TCH" + assert "https://stac-extensions.github.io/grid/v1.1.0/schema.json" in item.stac_extensions + + def test_builds_from_non_consolidated_store(tmp_path: Path) -> None: """Regression: the builder must read a store that lacks root consolidated metadata. diff --git a/tests/test_s1_stac_per_acquisition.py b/tests/test_s1_stac_per_acquisition.py index 8411e9eb..893b8675 100644 --- a/tests/test_s1_stac_per_acquisition.py +++ b/tests/test_s1_stac_per_acquisition.py @@ -233,6 +233,14 @@ def test_per_slice_platform_and_no_datacube(tmp_path: Path) -> None: assert "created" not in item.properties +def test_grid_code_inherited_from_cube(tmp_path: Path) -> None: + """Per-acquisition items inherit the cube's grid extension + grid:code (the MGRS tile).""" + store = _make_acq_cube(tmp_path, {"descending": [(T_EARLY, "S1A")]}) + item = build_s1_rtc_per_acquisition_items(store, orbit="descending", collection_id="acq")[0] + assert item.properties["grid:code"] == "MGRS-31TCH" + assert "https://stac-extensions.github.io/grid/v1.1.0/schema.json" in item.stac_extensions + + def test_footprint_is_run_orbit_not_cube_union(tmp_path: Path) -> None: """A per-acq item must carry ITS orbit's footprint, not the cube's union of both orbits' extents.""" # ascending shifted east so the cube union bbox is wider than the descending orbit alone. From 0eb157ae786442f20d94d26e3db95e2d0173effc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Houpert?= <10154151+lhoupert@users.noreply.github.com> Date: Tue, 7 Jul 2026 15:17:07 +0100 Subject: [PATCH 19/19] chore: regenerate uv.lock with pystac after rebase onto main The rebase took main's uv.lock (--ours) at the stac-builder conflict; regenerate so the lockfile includes pystac (required by the S1 RTC STAC builder) on top of main's zarr-cm 0.4.1 / pyright dependency bump. Co-Authored-By: Claude Fable 5 --- uv.lock | 2 ++ 1 file changed, 2 insertions(+) diff --git a/uv.lock b/uv.lock index a7c71ef4..d4500c9c 100644 --- a/uv.lock +++ b/uv.lock @@ -803,6 +803,7 @@ dependencies = [ { name = "pydantic" }, { name = "pydantic-zarr" }, { name = "pyproj" }, + { name = "pystac" }, { name = "rioxarray" }, { name = "s3fs" }, { name = "structlog" }, @@ -849,6 +850,7 @@ requires-dist = [ { name = "pydantic", specifier = ">=2.12" }, { name = "pydantic-zarr", specifier = ">=0.8.0" }, { name = "pyproj", specifier = ">=3.7.0" }, + { name = "pystac", specifier = ">=1.8.0" }, { name = "rioxarray", specifier = ">=0.13.0" }, { name = "s3fs", specifier = ">=2024.6.0" }, { name = "structlog", specifier = ">=25.5.0" },