diff --git a/.github/prompts/phase2-geotiff-ingestion-plan.md b/.github/prompts/phase2-geotiff-ingestion-plan.md new file mode 100644 index 00000000..ffd058e2 --- /dev/null +++ b/.github/prompts/phase2-geotiff-ingestion-plan.md @@ -0,0 +1,471 @@ +# Phase 2 — GeoTIFF Ingestion: Detailed Implementation Plan + +## Objective + +Productionise the prototype GeoTIFF → Zarr V3 ingestion pipeline (`analysis/s1_grd_rtc_prototype.py` and `analysis/s1_real_geotiff_to_zarr.py`) into a production module at `src/eopf_geozarr/conversion/s1_ingest.py` with proper logging, error handling, and test coverage. + +Phase 2 covers **data ingestion only** (not conditions, overviews, CLI, or S3) — those are Phases 3–4. However, the internal design should cleanly accommodate them. + +--- + +## Prior Art — What Exists + +| Asset | Path | Relevance | +|-------|------|-----------| +| Phase 0 synthetic prototype | `analysis/s1_grd_rtc_prototype.py` | ~900 lines. Creates store, appends 2 acquisitions with overviews, validates, consolidates metadata. Creates 1D spatial coordinate arrays (`x`, `y`) at every resolution level. All synthetic 256×256 data. | +| Phase 0 real-data script | `analysis/s1_real_geotiff_to_zarr.py` | 700 lines. Handles real S1Tiling filenames, datetime normalisation (`2025:02:10T06:09:20Z`), gamma_area conditions, validation report. Validated on 3 real 10980×10980 acquisitions. | +| Phase 1 data model | `src/eopf_geozarr/data_api/s1_rtc.py` | ~350 lines. Pydantic-zarr V3 models: `S1RtcRoot`, `S1RtcOrbitGroup`, `S1RtcNativeResolutionDataset`, etc. Strict validation via `@model_validator`. Now uses typed `Multiscales`/`MultiscalesScaleLevel` Pydantic models instead of `dict[str, Any]`. Validates `spatial:dimensions == ["y", "x"]` (lowercase). | +| Phase 1 test fixture | `tests/_test_data/s1_rtc_examples/s1-grd-rtc-31TCH.json` | Full JSON metadata for a 3-acquisition store with conditions. | +| Existing utilities | `src/eopf_geozarr/conversion/utils.py` | `calculate_aligned_chunk_size()`, `downsample_2d_array()` — **reuse directly** | +| Existing FS abstraction | `src/eopf_geozarr/conversion/fs_utils.py` | S3/local filesystem helpers — **reuse for S3 path support in Phase 4** | +| Real GeoTIFF metadata | `analysis/s1tiling_output_metadata.json` | Captures all 32 tags from actual S1Tiling 1.4.0 outputs | +| Zarr conventions library | `zarr_cm` (geo_proj, multiscales, spatial) | UUIDs and schema URLs used in Phase 1 models — **reuse for constants** | + +### Key decisions already validated in Phase 0 + +- `zarr.create_array(shards=, compressors=)` works (NOT `codecs=`) +- `array.resize()` works on sharded arrays, preserves existing data +- `rasterio src.transform` returns Affine ordering natively (no GDAL conversion needed) +- Inner chunk size must evenly divide shard: `calculate_aligned_chunk_size(10980, 512)` → 366 +- S1Tiling datetime format `"2025:02:10T06:09:20Z"` needs normalisation (colons in date part) +- border_mask is per-polarisation; prototype uses VV mask as primary +- Overview ceiling division: `ceil(10980/2) = 5490`, `ceil(5490/3) = 1830`, etc. +- **Lowercase dimension names**: `["y", "x"]` throughout (not `["Y", "X"]`) — required by titiler-eopf +- **1D spatial coordinate arrays** (`x`, `y`) must exist at every resolution level — required by GeoZarr readers +- **Consolidated metadata** at orbit and root levels after all ingestions — single-request metadata reads +- **`create_array(data=...)` cannot be combined with `dtype=`** — zarr-python 3.1.1 raises ValueError + +--- + +## Architecture + +### Module: `src/eopf_geozarr/conversion/s1_ingest.py` + +One file, ~500 lines estimated. Internal organisation: + +``` +s1_ingest.py +├── Constants (conventions, overview chain) +├── GeoTIFF metadata extraction +│ ├── S1TilingMetadata (dataclass) +│ ├── S1TILING_FILENAME_PATTERN (regex) +│ ├── extract_geotiff_metadata() +│ └── parse_s1tiling_filename() +├── Store creation +│ ├── compute_multiscales_layout() +│ ├── _create_spatial_coordinate_arrays() +│ └── create_s1_store() +├── Acquisition ingestion +│ ├── ingest_s1tiling_acquisition() ← PUBLIC API +│ └── _normalise_s1tiling_datetime() +├── Consolidation +│ └── consolidate_s1_store() ← PUBLIC API +├── File discovery +│ └── discover_s1tiling_acquisitions() ← PUBLIC API +└── (future hooks for conditions and overview levels) +``` + +### Tests: `tests/test_s1_rtc_ingest.py` + +~300 lines estimated. Uses synthetic GeoTIFFs created via `rasterio` in fixtures. + +--- + +## Step-by-step Implementation + +### Step 1: Constants and dataclass + +**What**: Define `S1TilingMetadata` dataclass and shared constants. + +**Details**: +- `S1TilingMetadata`: A frozen dataclass holding all fields extracted from a GeoTIFF (crs, spatial_transform, shape, bounds, datetime, absolute_orbit, relative_orbit, platform, calibration, input_s1_images). NOT a Pydantic model — this is simple data transfer, not validation. +- `OVERVIEW_CHAIN`: same `[("r10m", None, 1), ("r20m", "r10m", 2), ...]` tuple list +- Zarr convention constants: import from `zarr_cm` (already used in Phase 1 models) rather than hardcoding UUIDs. Specifically: + ```python + from zarr_cm import geo_proj, multiscales as multiscales_cm, spatial as spatial_cm + MULTISCALES_UUID = multiscales_cm.UUID + GEO_PROJ_UUID = geo_proj.UUID + SPATIAL_UUID = spatial_cm.UUID + ``` +- `ZARR_CONVENTIONS` list: build from `zarr_cm` library attributes (schema_url, spec_url, name, description, uuid) rather than hardcoding the full dicts. Check if `zarr_cm` exposes these as structured objects; if not, hardcode as in prototype. +- `S1TILING_FILENAME_PATTERN`: the regex from `s1_real_geotiff_to_zarr.py` + +**Reuse**: Constants pattern from `s1_rtc.py` (Phase 1). `calculate_aligned_chunk_size` from `utils.py`. + +--- + +### Step 2: Metadata extraction + +**What**: `extract_geotiff_metadata(path) -> S1TilingMetadata` and `parse_s1tiling_filename(filename) -> dict`. + +**Details**: + +`extract_geotiff_metadata`: +- Opens GeoTIFF with `rasterio.open()` (read-only) +- Reads CRS (`str(src.crs)`), transform (Affine → list of 6 floats `[a, b, c, d, e, f]`), bounds, shape +- Reads custom tags: `ACQUISITION_DATETIME`, `ORBIT_NUMBER`, `RELATIVE_ORBIT_NUMBER`, `FLYING_UNIT_CODE`, `CALIBRATION`, `INPUT_S1_IMAGES` +- Normalises datetime: `_normalise_s1tiling_datetime("2025:02:10T06:09:20Z")` → `"2025-02-10T06:09:20"` +- Returns `S1TilingMetadata` instance +- Raises `ValueError` if critical tags are missing (ACQUISITION_DATETIME, ORBIT_NUMBER, RELATIVE_ORBIT_NUMBER, FLYING_UNIT_CODE) +- structlog for info-level logging of extracted metadata + +`parse_s1tiling_filename`: +- Applies `S1TILING_FILENAME_PATTERN` regex to filename string +- Returns dict with keys: platform, tile, pol, orbit_dir, rel_orbit, acq_stamp, is_mask +- Returns `None` if filename doesn't match (not an error — allows callers to skip) + +`_normalise_s1tiling_datetime`: +- Input: `"2025:02:10T06:09:20Z"` (S1Tiling format with colons in date) +- Output: `"2025-02-10T06:09:20"` (ISO 8601 for `np.datetime64`) +- Strip trailing Z, split on T, replace colons with hyphens in date part only + +**Source to lift from**: `s1_real_geotiff_to_zarr.py` lines 107–132 (extract_geotiff_metadata) and lines 91–105 (S1TILING_PATTERN + parse logic) and lines 392–399 (datetime normalisation). + +--- + +### Step 3: Multiscales layout computation + +**What**: `compute_multiscales_layout(native_shape, native_transform) -> list[dict]` + +**Details**: +- Takes native shape `[H, W]` and transform `[a, b, c, d, e, f]` (Affine ordering) +- Iterates OVERVIEW_CHAIN, producing one layout entry per level +- Each entry: `{"asset": name, "spatial:shape": [h, w], "spatial:transform": [a,b,c,d,e,f], "transform": {"scale": [f, f]}, ...}` +- Native level has `"transform": {"scale": [1.0, 1.0]}`; deeper levels have `"derived_from": parent_name` +- Shape computation: `ceil(parent_h / factor)`, `ceil(parent_w / factor)` +- Transform: scale pixel size by factor (`transform[0] *= factor`, `transform[4] *= factor`), keep origin + +**Source to lift from**: `s1_grd_rtc_prototype.py` lines 129–173 (`compute_multiscales_layout`). Unchanged logic, production-ready. + +--- + +### Step 4: Store creation + +**What**: `create_s1_store(store_path, orbit_direction, metadata) -> zarr.Group` + +**Details**: +- Creates Zarr V3 store at `store_path` (local path string or Path) +- `zarr.open_group(str(store_path), mode="w", zarr_format=3)` +- Creates orbit direction group with full conventions attributes: + - `zarr_conventions`: ZARR_CONVENTIONS list + - `multiscales`: `{"layout": [...], "resampling_method": "average"}` + - `proj:code`: from metadata CRS + - `spatial:dimensions`: `["y", "x"]` (**lowercase** — required by titiler-eopf) + - `spatial:bbox`: from metadata bounds +- For each level in layout: + - Create level group with `spatial:shape` and `spatial:transform` attributes + - Create `vv`, `vh`, `border_mask` arrays with: + - `shape=(0, level_h, level_w)` — time axis starts at 0 + - `dtype`: float32 for vv/vh, uint8 for border_mask + - `fill_value`: NaN for float32, 0 for uint8 + - `chunks`: `(1, best_chunk(level_h), best_chunk(level_w))` using `calculate_aligned_chunk_size()` + - `shards`: `(1, level_h, level_w)` — one shard per timestep + - `compressors`: `BloscCodec(cname="zstd", clevel=5)` + - `dimension_names`: `["time", "y", "x"]` (**lowercase**) + - Create **1D spatial coordinate arrays** (`x`, `y`) at this level: + - Compute from the level's `spatial:transform` `[a, b, c, d, e, f]`: + - `x_coords = np.linspace(x_origin, x_origin + level_w * pixel_w, level_w, endpoint=False)` + - `y_coords = np.linspace(y_origin, y_origin + level_h * pixel_h, level_h, endpoint=False)` + - Create with `data=` param (not `shape=` + write), omit `dtype=` (inferred from data) + - `dimension_names`: `["x"]` and `["y"]` respectively + - Attributes: `units: "m"`, `long_name`, `standard_name: "projection_x_coordinate"` / `"projection_y_coordinate"`, `_ARRAY_DIMENSIONS: ["x"]` / `["y"]` + - **API note**: `create_array(data=...)` cannot specify `dtype=` — it is inferred +- Create coordinate variables at r10m only: + - `time`: int64, shape (0,), chunks (512,), dim_names ["time"] + - `absolute_orbit`: int32, same + - `relative_orbit`: int32, same + - `platform`: ` int` + +This is the **main public API** for Phase 2. + +**Details**: + +1. **Extract metadata** from vv_path via `extract_geotiff_metadata()` +2. **Create-or-open store**: + - If store doesn't exist → call `create_s1_store()` + - If store exists → open in `mode="r+"` + - If store exists but orbit direction group doesn't exist → create it +3. **Read GeoTIFF pixel data**: + - `rasterio.open(vv_path).read(1)` → vv_data (float32) + - `rasterio.open(vh_path).read(1)` → vh_data (float32) + - `rasterio.open(border_mask_path).read(1).astype(np.uint8)` → mask_data +4. **Determine time index**: `current_size = r10m["vv"].shape[0]`, `new_size = current_size + 1` +5. **Generate overviews** from native data: + - For each level in `OVERVIEW_CHAIN[1:]`: downsample from previous level + - vv/vh: `downsample_2d` with `"average"` method + - border_mask: `downsample_2d` with `"nearest"` method + - Use existing `downsample_2d_array()` from `utils.py` (different signature — takes target_height/target_width instead of factor, so compute target sizes first using ceiling division) + - **NOTE**: The existing `downsample_2d_array()` in `utils.py` has a subtly different interface from the prototype's `downsample_2d()`. The prototype takes a `factor` and computes target dims internally with ceiling division. The production util takes target dims directly. Need to compute `target_h = ceil(h / factor)`, `target_w = ceil(w / factor)` before calling. + - **ALTERNATIVE**: The prototype's `downsample_2d()` is simpler and handles edge padding correctly for non-divisible sizes. Consider adding it as a private helper in `s1_ingest.py` or adding a `factor`-based variant to `utils.py`. Decision: **Use the prototype's `downsample_2d()` as a private helper** — it's purpose-built for power-of-2/3 downsampling with padding, which is the S1 use case. The existing `downsample_2d_array()` was designed for arbitrary target dimensions (S2 overview path). +6. **Write data at all levels**: + - For each `(level_name, data_tuple)`: + - `level["vv"].resize((new_size, h, w))` + - `level["vv"][current_size, :, :] = vv_lev` + - Same for vh, border_mask +7. **Append coordinate variables**: + - Resize all 1-D arrays to `new_size` + - Write `time[current_size] = np.datetime64(dt).astype("datetime64[ns]").astype(np.int64)` + - Write `absolute_orbit[current_size] = meta.absolute_orbit` + - Write `relative_orbit[current_size] = meta.relative_orbit` + - Write `platform[current_size] = meta.platform` +8. **Return** the time index (`current_size`) + +**Logging** (structlog): +- INFO: "Ingesting S1 acquisition", vv_path, orbit_direction, time_index +- INFO: "GeoTIFF read complete", read_time_s, vv_min, vv_max, mask_coverage_pct +- INFO: "Overviews generated", overview_time_s +- INFO: "Zarr write complete", write_time_s, levels_written + +**Error handling** (at system boundary — GeoTIFF files are external input): +- `FileNotFoundError`: if vv/vh/mask paths don't exist +- `ValueError`: if GeoTIFF CRS doesn't match store's `proj:code` (on append) +- `ValueError`: if GeoTIFF shape doesn't match store's native shape (on append) +- Let zarr exceptions propagate naturally for I/O errors + +**Source to lift from**: `s1_real_geotiff_to_zarr.py` lines 352–432 (`ingest_acquisition`). Add create-or-open logic, CRS/shape consistency checks, structlog. + +--- + +### Step 6: File discovery utility + +**What**: `discover_s1tiling_acquisitions(input_dir) -> list[dict]` + +**Details**: +- Glob `input_dir/*.tif`, apply `parse_s1tiling_filename()` to each +- Group by `(platform, tile, orbit_dir, rel_orbit, acq_stamp)` key +- Validate each group has vv, vh, vv_mask, vh_mask +- Return list of dicts with keys: `platform, tile, orbit_dir, rel_orbit, acq_stamp, vv, vh, vv_mask, vh_mask` (paths as `Path` objects) +- Log warnings for incomplete acquisitions + +This is mainly useful for the CLI batch ingest (Phase 4) but is simple to implement and test now. + +**Source to lift from**: `s1_real_geotiff_to_zarr.py` lines 142–178 (`discover_acquisitions`). + +--- + +### Step 7: Module exports and wiring + +**What**: Register public API in `__init__.py`. + +**Details**: +- Add to `src/eopf_geozarr/conversion/__init__.py`: + ```python + from .s1_ingest import ( + ingest_s1tiling_acquisition, + discover_s1tiling_acquisitions, + extract_geotiff_metadata, + ) + ``` +- Add to `__all__` + +--- + +### Step 8: Test — synthetic GeoTIFF fixtures + +**What**: Create pytest fixtures that produce synthetic S1Tiling GeoTIFFs. + +**File**: `tests/test_s1_rtc_ingest.py` + +**Fixtures**: +- `s1_geotiff_dir(tmp_path)`: Creates a temp directory with synthetic GeoTIFFs for 2 acquisitions: + - `s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC.tif` + - `s1a_32TQM_vh_ASC_037_20230115t061234_GammaNaughtRTC.tif` + - `s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC_BorderMask.tif` + - `s1a_32TQM_vh_ASC_037_20230115t061234_GammaNaughtRTC_BorderMask.tif` + - (second acquisition: `20230127t061235`) + - All 256×256, EPSG:32633, with proper tags (`ACQUISITION_DATETIME`, `ORBIT_NUMBER`, etc.) +- `s1_store_path(tmp_path)`: Returns a clean path for Zarr store output + +**Helper**: `_create_synthetic_geotiff(path, data, crs, transform, tags)` — same as prototype's `create_test_geotiff()`. + +--- + +### Step 9: Test — metadata extraction + +**Tests**: +- `test_extract_geotiff_metadata()`: Verify all fields populated from tags +- `test_extract_geotiff_metadata_normalises_datetime()`: `"2025:02:10T06:09:20Z"` → `"2025-02-10T06:09:20"` +- `test_extract_geotiff_metadata_raises_on_missing_tags()`: Missing ORBIT_NUMBER → ValueError +- `test_parse_s1tiling_filename()`: Correct field extraction for VV, VH, mask variants +- `test_parse_s1tiling_filename_returns_none_for_unknown()`: Non-matching filename + +--- + +### Step 10: Test — store creation + +**Tests**: +- `test_create_s1_store_structure()`: Creates store, verifies group hierarchy (root → ascending → r10m…r720m) +- `test_create_s1_store_conventions()`: Verifies `zarr_conventions`, `proj:code`, `spatial:dimensions == ["y", "x"]`, `spatial:bbox` on orbit group +- `test_create_s1_store_array_metadata()`: Verifies `dimension_names == ["time", "y", "x"]`, dtype, fill_value, shape=(0,H,W) for data arrays +- `test_create_s1_store_coordinate_arrays()`: Verifies time, absolute_orbit, relative_orbit, platform at r10m only +- `test_create_s1_store_overview_shapes()`: Verifies consistent shape chain (ceiling division) +- `test_create_s1_store_spatial_coordinate_arrays()`: Verifies `x` and `y` 1D arrays exist at EVERY resolution level, with correct shape, attributes (`units`, `standard_name`, `_ARRAY_DIMENSIONS`), and values derived from spatial:transform + +--- + +### Step 11: Test — ingestion (create + append) + +**Tests**: +- `test_ingest_first_acquisition()`: Ingest into non-existing store → store created, time_index=0, data readable +- `test_ingest_second_acquisition_appends()`: Ingest twice → shape[0]=2, both timesteps have data +- `test_ingest_preserves_data_integrity()`: Write known data, read back, compare values (within float32 tolerance for vv/vh, exact for mask) +- `test_ingest_coordinate_values()`: Verify time, orbit, platform values written correctly +- `test_ingest_overview_consistency()`: Verify overview shapes follow ceiling division chain +- `test_ingest_rejects_mismatched_crs()`: Second acquisition with different CRS → ValueError +- `test_ingest_rejects_mismatched_shape()`: Second acquisition with different shape → ValueError +- `test_ingest_xarray_roundtrip()`: After 2 ingestions, `xr.open_zarr(r10m_path)` succeeds, `ds.sortby("time")` works + +--- + +### Step 12b: Test — consolidation + +**Tests**: +- `test_consolidate_s1_store()`: After ingestion + `consolidate_s1_store()`, orbit group and root have `consolidated_metadata` present +- `test_consolidate_after_all_ingestions()`: Verify that consolidation after 2 ingestions records the correct array shapes (not stale from mid-ingestion) + +--- + +### Step 12: Test — file discovery + +**Tests**: +- `test_discover_acquisitions_groups_correctly()`: 2 acquisitions × 4 files each → 2 groups +- `test_discover_acquisitions_warns_on_incomplete()`: Missing VH file → warning logged, still returns partial +- `test_discover_acquisitions_skips_non_matching()`: Random .tif files in directory → ignored + +--- + +## What is NOT in Phase 2 + +These are explicitly deferred: + +| Item | Phase | Reason | +|------|-------|--------| +| Conditions ingestion (gamma_area, LIA) | Phase 3 | Separate data path, no time dimension | +| Conditions group with own conventions | Phase 3 | Depends on conditions ingestion | +| CLI subcommands (ingest-s1, etc.) | Phase 4 | Needs stable public API first | +| S3 output support | Phase 4 | Requires fs_utils integration; local-first | +| Multiscale generation as separate concern | Phase 3 | Currently embedded in ingest; may refactor | +| STAC item creation/update | Phase 5 | Needs STAC collection definition first | +| Validation via Phase 1 Pydantic models | Phase 3+ | Optional; store structure is validated by tests | +| border_mask combining (VV ∪ VH) | Future | Prototype uses VV-only; acceptable for now | +| Phase 1 model update for x/y members | Phase 2 Follow-up | Add `x`, `y` to `S1RtcOverviewResolutionMembers` and `S1RtcNativeResolutionMembers` TypedDicts | + +--- + +## Dependency/Import Map + +``` +s1_ingest.py imports: + ├── from __future__ import annotations + ├── dataclasses (dataclass, frozen) + ├── math (ceil) + ├── pathlib (Path) + ├── re + ├── numpy as np + ├── rasterio + ├── structlog + ├── zarr + ├── zarr.codecs (BloscCodec) + ├── zarr_cm (geo_proj, multiscales, spatial) # UUIDs and schema URLs + └── eopf_geozarr.conversion.utils (calculate_aligned_chunk_size) + +test_s1_rtc_ingest.py imports: + ├── xarray as xr # for roundtrip tests + └── (all of the above transitively via the module under test) +``` + +No new external dependencies. All imports already in the project. + +--- + +## Estimated Scope + +| Component | Lines (est.) | Complexity | +|-----------|-------------|------------| +| Constants + S1TilingMetadata | ~50 | Low | +| extract_geotiff_metadata + helpers | ~70 | Low | +| compute_multiscales_layout | ~50 | Low (direct lift) | +| _create_spatial_coordinate_arrays | ~40 | Low (direct lift from prototype) | +| create_s1_store | ~100 | Medium (now includes x/y coord creation) | +| ingest_s1tiling_acquisition | ~120 | Medium-High | +| consolidate_s1_store | ~15 | Low | +| discover_s1tiling_acquisitions | ~50 | Low | +| _downsample_2d (private helper) | ~30 | Low (direct lift) | +| **Total s1_ingest.py** | **~525** | | +| Test fixtures | ~80 | Low | +| Test cases (15 tests) | ~300 | Medium | +| **Total test_s1_rtc_ingest.py** | **~380** | | + +--- + +## Implementation Order for the Executing Agent + +The steps above are ordered for incremental, testable progress. A reasonable execution sequence: + +1. **Steps 1–3**: Constants, metadata extraction, multiscales layout — pure functions, no I/O beyond rasterio reads +2. **Step 8**: Create test fixtures (needed for all subsequent testing) +3. **Steps 9**: Test metadata extraction +4. **Step 4**: Store creation +5. **Step 10**: Test store creation +6. **Step 5**: Acquisition ingestion +7. **Steps 6**: File discovery +8. **Steps 11–12**: Integration tests +9. **Step 7**: Wire up exports + +Each step can be committed independently. The agent should run `pytest tests/test_s1_rtc_ingest.py -v` after each test step to verify green. +--- + +## Appendix A: Post-Phase-2 Model Fix + +The Phase 1 `S1RtcOverviewResolutionMembers` and `S1RtcNativeResolutionMembers` TypedDicts do not declare `x` and `y` members, yet the store structure requires 1D spatial coordinate arrays at every level. After Phase 2 is complete and tested, submit a follow-up PR to `s1_rtc.py` adding: + +```python +class S1RtcNativeResolutionMembers(TypedDict, closed=True, total=False): + # existing: vv, vh, border_mask, time, absolute_orbit, relative_orbit, platform + x: ArraySpec[Any] # (x,) float64 + y: ArraySpec[Any] # (y,) float64 + +class S1RtcOverviewResolutionMembers(TypedDict, closed=True, total=False): + # existing: vv, vh, border_mask + x: ArraySpec[Any] # (x,) float64 + y: ArraySpec[Any] # (y,) float64 +``` + +Update the test fixture JSON (`s1-grd-rtc-31TCH.json`) and tests to include `x`/`y` in the members. + +--- + +## Appendix B: Consolidation Step + +Consolidation is NOT called inside `ingest_s1tiling_acquisition()` — it must be called **after all ingestions for a batch** via the separate `consolidate_s1_store()` function. This is because consolidated metadata caches array shapes; consolidating mid-flow causes `BoundsCheckError` on subsequent `resize()` calls. + +```python +def consolidate_s1_store(store_path: str | Path, orbit_direction: str) -> None: + """Consolidate metadata at orbit direction and root levels. + + Must be called AFTER all ingestions complete. + """ + zarr.consolidate_metadata(str(store_path), path=orbit_direction, zarr_format=3) + zarr.consolidate_metadata(str(store_path), zarr_format=3) +``` + +--- + +## Appendix C: Dimension Name Convention + +All dimension names and `spatial:dimensions` MUST use **lowercase** `["y", "x"]` and `["time", "y", "x"]`, not uppercase. This was discovered during post-validation testing with titiler-eopf and is now the canonical convention across the prototype, Phase 1 models, and implementation plan. + +**Note**: The `s1_real_geotiff_to_zarr.py` analysis script still has some uppercase `["Y", "X"]` references (in `create_s1_store` and `ingest_gamma_area`). These are inconsistencies in the analysis script only; the production code must use lowercase. \ No newline at end of file diff --git a/.github/prompts/s1-grd-rtc-implementation-plan-v2.md b/.github/prompts/s1-grd-rtc-implementation-plan-v2.md new file mode 100644 index 00000000..7133fbeb --- /dev/null +++ b/.github/prompts/s1-grd-rtc-implementation-plan-v2.md @@ -0,0 +1,823 @@ +# Sentinel-1 GRD γ0T RTC — Implementation Plan for data-model codebase + +## Context + +EOPF Explorer needs Sentinel-1 GRD data in its catalog. ESA has not provided NRB ARCO products. We use **S1Tiling** (CNES/OTB) to produce γ0T RTC GeoTIFFs on the Sentinel-2 MGRS grid, then convert those GeoTIFFs into GeoZarr stores following the EOPF hierarchy. + +This plan is scoped to the `EOPF-Explorer/data-model` repository. It extends the existing codebase (which handles S2 EOPF Zarr → GeoZarr conversion) with a new ingestion path for S1Tiling GeoTIFF outputs. + +Repository: https://github.com/EOPF-Explorer/data-model +Current version: v0.8.0 +GeoZarr spec: **GeoZarr 1.0** — modular Zarr Conventions framework (post-Rome Summit Dec 2025) +Mini-spec: https://eopf-explorer.github.io/data-model/geozarr-minispec/ +Key dependency: S1Tiling 1.4.0 (https://s1-tiling.pages.orfeo-toolbox.org/s1tiling/latest/) + +### Architecture decision: two-step pipeline, GeoTIFF as handoff + +S1Tiling is treated as an **external black box** that produces GeoTIFFs. It is NOT +integrated as a Python library dependency in this repository. + +Rationale: S1Tiling orchestrates OTB C++ applications that fundamentally produce +GeoTIFF on disk (SAFE → OTB → GeoTIFF). Even calling `s1_process()` from Python, +the pixel data still materialises as GeoTIFF intermediates. Embedding OTB (a ~1GB +C++ binary with its own GDAL, proj, and LD_LIBRARY_PATH) into our image would add +dependency risk with no performance benefit — the GeoTIFF intermediate exists +regardless. + +The pipeline is therefore two Argo workflow steps: + +``` +Step 1: cnes/s1tiling Docker → GeoTIFFs on shared volume / S3 +Step 2: eopf-geozarr Docker → reads GeoTIFFs, writes Zarr + STAC +``` + +This repository owns Step 2 only. S1Tiling configuration and Docker image +management are pipeline concerns (Argo workflow YAML), not data-model concerns. + +--- + +## 1. GeoZarr Conventions Used + +Three core Zarr conventions declared in `zarr_conventions` at each relevant group level: + +| Convention | Namespace | UUID | Purpose | +|------------|-----------|------|---------| +| multiscales | `multiscales` | `d35379db-88df-4056-af3a-620245f8e347` | Pyramid layout | +| geo-proj | `proj:` | `f17cb550-5864-4468-aeb7-f3180cfb622f` | CRS encoding | +| spatial | `spatial:` | `689b58e2-cf7b-45e0-9fff-9cfc0883d6b4` | Array index → spatial coordinate | + +Key principles from the mini-spec: +- **No `grid_mapping` 0D arrays** — CRS via `proj:code` at group level, inherited by child arrays +- **No mandatory CF conventions** — CF metadata is compatible but not required +- **`dimension_names`** in Zarr V3 array metadata (not `_ARRAY_DIMENSIONS` attribute) +- **`spatial:transform`** in Rasterio/Affine ordering `[a, b, c, d, e, f]` (NOT GDAL ordering) +- **Multiscales** use `layout` array with `asset`, `derived_from`, `transform.scale` +- **1D coordinate arrays** (`x`, `y`) must exist at every resolution level for GeoZarr reader compatibility (e.g. titiler-eopf) +- **Consolidated metadata** at root and orbit direction group levels for performant reads + +--- + +## 2. Zarr Store Structure + +One Zarr V3 store per S2 MGRS tile. Store name: `s1-grd-rtc-{tile_id}.zarr` + +``` +s1-grd-rtc-32TQM.zarr/ +├── zarr.json # Root group (no conventions at root) +│ +├── ascending/ +│ ├── zarr.json # Group: zarr_conventions [multiscales, proj:, spatial:] +│ │ # proj:code: "EPSG:32633" +│ │ # spatial:dimensions: ["y", "x"] +│ │ # spatial:bbox: [xmin, ymin, xmax, ymax] +│ │ # multiscales: { layout: [...] } +│ │ +│ ├── r10m/ # Native resolution dataset (asset: "r10m") +│ │ ├── zarr.json # Group: spatial:shape, spatial:transform +│ │ ├── vv/ # (time, y, x) float32 +│ │ │ ├── zarr.json # dimension_names: ["time", "y", "x"] +│ │ │ └── c/{t}/0/0 # One shard per timestep +│ │ ├── vh/ # (time, y, x) float32 +│ │ │ └── c/{t}/0/0 +│ │ ├── border_mask/ # (time, y, x) uint8 +│ │ │ └── c/{t}/0/0 +│ │ ├── x/ # (x,) float64 — projection x coordinates +│ │ ├── y/ # (y,) float64 — projection y coordinates +│ │ ├── time/ # (time,) datetime64[ns] +│ │ ├── absolute_orbit/ # (time,) int32 +│ │ ├── relative_orbit/ # (time,) int32 +│ │ └── platform/ # (time,) str +│ │ +│ ├── r20m/ # Overview level 1 (2x from r10m) +│ │ ├── zarr.json # spatial:shape: [5490, 5490] +│ │ │ # spatial:transform: [20.0, 0.0, ...] +│ │ ├── vv/ +│ │ ├── vh/ +│ │ ├── border_mask/ +│ │ ├── x/ # (x,) float64 — projection x at this resolution +│ │ └── y/ # (y,) float64 — projection y at this resolution +│ │ +│ ├── r60m/ # Overview level 2 (3x from r20m) +│ ├── r120m/ # Overview level 3 (2x from r60m) +│ ├── r360m/ # Overview level 4 (3x from r120m) +│ ├── r720m/ # Overview level 5 (2x from r360m) +│ │ +│ └── conditions/ # Time-invariant, NOT in multiscales layout +│ ├── zarr.json # proj:code, spatial:dimensions, spatial:transform +│ ├── gamma_area_{orbit}/ # (Y, X) float32 — one per relative orbit number +│ ├── lia_{orbit}/ # (Y, X) float32 — sin(LIA) [if available] +│ └── incidence_angle_{orbit}/ # (Y, X) float32 [if available] +│ +└── descending/ + └── (same structure) +``` + +### Key design decisions + +**Multiscale levels are named by resolution** (`r10m`, `r20m`, `r60m`, ...) to match the S2 convention used in the existing codebase. Each level is a Dataset containing the same set of variables. + +**Coordinate variables** (`time`, `absolute_orbit`, `relative_orbit`, `platform`) live inside the native resolution dataset (`r10m/`) because they follow the Dataset rule: for each dimension name in a data variable, there must be a matching 1D coordinate variable. Overview levels share the same time dimension but don't need separate coordinate copies (they reference the same time axis). + +**1D spatial coordinate arrays** (`x`, `y`) are required at **every resolution level** (r10m through r720m). They are computed from the level's `spatial:transform` using `np.linspace` and carry CF-standard attributes (`units: "m"`, `standard_name: "projection_x_coordinate"` / `"projection_y_coordinate"`). Without these arrays, GeoZarr readers like titiler-eopf cannot resolve spatial coordinates from the data. + +**Consolidated metadata** is written at the root and orbit direction group levels after all ingestions are complete. This embeds the full hierarchy metadata in the group-level `zarr.json`, enabling single-request metadata reads. Important: consolidation must happen **after** all resize/append operations — consolidating mid-ingestion caches stale array shapes and breaks subsequent writes. + +**Conditions** sit outside the multiscales layout as a separate group. They are (Y, X) only — no time dimension — and are per orbit, not per acquisition. They carry their own `proj:` and `spatial:` conventions. Each array is named with the relative orbit number suffix (e.g. `gamma_area_008`). Only `gamma_area` is confirmed as a direct S1Tiling output; `lia` and `incidence_angle` may require additional S1Tiling configuration or post-processing to produce as separate files. + +**border_mask** is included as a variable alongside vv/vh in each resolution level. It shares the (time, Y, X) shape and gets downsampled with the overviews (using `nearest` resampling for masks, not `average`). + +### Chunking and sharding + +- **Chunks**: `(1, C, C)` where C = largest divisor of tile dimension ≤ 512 (e.g. 366 for 10980, since 10980/30=366). Zarr sharding requires inner chunks to **evenly divide** the shard; 512 does NOT divide 10980. +- **Shards**: `(1, H, W)` — one shard file = one timestep = all spatial chunks +- **Physical files**: `vv/c/{time_index}/0/0` — each new acquisition → one new shard file per array +- **chunk_key_encoding**: `{"name": "default", "configuration": {"separator": "/"}}` +- **Overview shards**: same strategy, naturally smaller per level + +### Time dimension — append model + +- Append-order integer index as dimension axis +- Actual datetime stored in `time` coordinate variable +- Non-monotonic time is expected and by design +- Consumers use `ds.sortby('time')` — no performance penalty since each timestep is an independent shard + +### Multiscales metadata (ascending/zarr.json) + +```json +{ + "zarr_conventions": [ + {"uuid": "d35379db-...", "name": "multiscales", ...}, + {"uuid": "f17cb550-...", "name": "proj:", ...}, + {"uuid": "689b58e2-...", "name": "spatial:", ...} + ], + "multiscales": { + "layout": [ + { + "asset": "r10m", + "transform": {"scale": [1.0, 1.0]}, + "spatial:shape": [10980, 10980], + "spatial:transform": [10.0, 0.0, 500000.0, 0.0, -10.0, 5000000.0] + }, + { + "asset": "r20m", + "derived_from": "r10m", + "transform": {"scale": [2.0, 2.0], "translation": [0.0, 0.0]}, + "spatial:shape": [5490, 5490], + "spatial:transform": [20.0, 0.0, 500000.0, 0.0, -20.0, 5000000.0] + }, + { + "asset": "r60m", + "derived_from": "r20m", + "transform": {"scale": [3.0, 3.0], "translation": [0.0, 0.0]}, + "spatial:shape": [1830, 1830], + "spatial:transform": [60.0, 0.0, 500000.0, 0.0, -60.0, 5000000.0] + }, + { + "asset": "r120m", + "derived_from": "r60m", + "transform": {"scale": [2.0, 2.0], "translation": [0.0, 0.0]}, + "spatial:shape": [915, 915], + "spatial:transform": [120.0, 0.0, 500000.0, 0.0, -120.0, 5000000.0] + }, + { + "asset": "r360m", + "derived_from": "r120m", + "transform": {"scale": [3.0, 3.0], "translation": [0.0, 0.0]}, + "spatial:shape": [305, 305], + "spatial:transform": [360.0, 0.0, 500000.0, 0.0, -360.0, 5000000.0] + }, + { + "asset": "r720m", + "derived_from": "r360m", + "transform": {"scale": [2.0, 2.0], "translation": [0.0, 0.0]}, + "spatial:shape": [153, 153], + "spatial:transform": [720.0, 0.0, 500000.0, 0.0, -720.0, 5000000.0] + } + ], + "resampling_method": "average" + }, + "proj:code": "EPSG:32633", + "spatial:dimensions": ["y", "x"], + "spatial:bbox": [500000.0, 4890200.0, 609800.0, 5000000.0] +} +``` + +### Array metadata example (r10m/vv/zarr.json) + +This is the actual on-disk format produced by zarr-python 3.1.1. Sharding is +represented as a codec wrapping the inner codecs. The `chunk_grid` holds the +**shard shape** (one shard = one timestep of the full spatial extent). The +**inner chunk shape** (`[1, 366, 366]`) lives inside the sharding codec config. + +```json +{ + "zarr_format": 3, + "node_type": "array", + "shape": [0, 10980, 10980], + "data_type": "float32", + "chunk_grid": { + "name": "regular", + "configuration": {"chunk_shape": [1, 10980, 10980]} + }, + "chunk_key_encoding": { + "name": "default", + "configuration": {"separator": "/"} + }, + "codecs": [ + { + "name": "sharding_indexed", + "configuration": { + "chunk_shape": [1, 366, 366], + "codecs": [ + {"name": "bytes", "configuration": {"endian": "little"}}, + {"name": "blosc", "configuration": {"cname": "zstd", "clevel": 5}} + ], + "index_codecs": [ + {"name": "bytes", "configuration": {"endian": "little"}}, + {"name": "crc32c"} + ], + "index_location": "end" + } + } + ], + "dimension_names": ["time", "y", "x"], + "fill_value": "NaN" +} +``` + +**Python API** (zarr-python 3.1.1): +```python +group.create_array( + "vv", shape=(0, 10980, 10980), dtype="float32", + chunks=(1, 366, 366), # inner chunk shape + shards=(1, 10980, 10980), # shard shape + compressors=zarr.codecs.BloscCodec(cname="zstd", clevel=5), + fill_value=float("nan"), + dimension_names=["time", "y", "x"], +) +``` + +Note: `shape[0]` starts at 0 and grows with each append. `dimension_names` is Zarr V3 native — no `_ARRAY_DIMENSIONS` attribute needed. Inner chunk size 366 = `best_chunk_size(10980)` = largest divisor of 10980 ≤ 512. + +--- + +## 3. New Modules + +### 3.1 Pydantic Models — `src/eopf_geozarr/models/sentinel1.py` + +This is the 404 page at https://eopf-explorer.github.io/data-model/models/sentinel1.md — it needs to be implemented. + +Extend existing base classes (inspect `models/sentinel2.py` for the pattern). The models should validate: + +- Store structure (ascending/descending groups) +- Zarr conventions declarations at each group level +- Array dimension_names consistency +- Coordinate variable existence for each dimension name in data variables +- Multiscales layout structure matching the mini-spec +- `proj:code`, `spatial:dimensions`, `spatial:transform` presence + +```python +# Core model sketch — refine against existing sentinel2.py base classes + +class S1GrdMeasurementsDataset(BaseModel): + """A single resolution level dataset (r10m, r20m, etc.).""" + vv: ArraySpec # (time, Y, X) float32 + vh: ArraySpec # (time, Y, X) float32 + border_mask: ArraySpec # (time, Y, X) uint8 + # Coordinate variables (at native resolution only) + time: Optional[ArraySpec] = None # (time,) datetime64 + absolute_orbit: Optional[ArraySpec] = None # (time,) int32 + relative_orbit: Optional[ArraySpec] = None # (time,) int32 + platform: Optional[ArraySpec] = None # (time,) str + +class S1GrdConditionsDataset(BaseModel): + """Time-invariant conditions per orbit. Arrays named with orbit suffix (e.g. gamma_area_008).""" + gamma_area: Dict[str, ArraySpec] # gamma_area_{orbit}: (Y, X) float32 — always present + lia: Optional[Dict[str, ArraySpec]] = None # lia_{orbit}: (Y, X) float32 — if available + incidence_angle: Optional[Dict[str, ArraySpec]] = None # incidence_angle_{orbit}: (Y, X) float32 — if available + +class S1GrdOrbitGroup(BaseModel): + """One orbit direction — carries multiscales, proj:, spatial: conventions.""" + multiscales: MultiscalesMetadata + proj_code: str # e.g., "EPSG:32633" + spatial_dimensions: List[str] # ["y", "x"] + spatial_bbox: List[float] + conditions: S1GrdConditionsDataset + +class S1GrdStore(BaseModel): + """Root store for one S2 MGRS tile.""" + tile_id: str + ascending: Optional[S1GrdOrbitGroup] = None + descending: Optional[S1GrdOrbitGroup] = None +``` + +### 3.2 GeoTIFF Ingestion — `src/eopf_geozarr/conversion/s1_ingest.py` + +**Public API:** + +```python +def ingest_s1tiling_acquisition( + vv_path: str, + vh_path: str, + mask_path: str, + zarr_store: str, + tile_id: str, + orbit_direction: str, # "ascending" or "descending" +) -> int: + """ + Append one S1Tiling acquisition to the Zarr store. + Creates the store with full zarr_conventions metadata if first acquisition. + Returns the time index of the appended acquisition. + + Metadata extraction from GeoTIFF tags: + - ACQUISITION_DATETIME → time coordinate + - ORBIT_NUMBER → absolute_orbit coordinate + - RELATIVE_ORBIT_NUMBER → relative_orbit coordinate + - FLYING_UNIT_CODE → platform coordinate + - CRS + GeoTransform → proj:code + spatial:transform (Rasterio ordering) + """ + +def ingest_s1tiling_conditions( + lia_path: str, + incidence_angle_path: str, + gamma_area_path: str, + zarr_store: str, + tile_id: str, + orbit_direction: str, +) -> None: + """ + Write time-invariant condition arrays. + Conditions group gets its own proj: and spatial: conventions. + """ +``` + +**Store creation** (first acquisition): +1. Write root `zarr.json` (minimal, no conventions at root) +2. Write `ascending/zarr.json` with full `zarr_conventions` array, `multiscales` layout, `proj:code`, `spatial:dimensions`, `spatial:bbox` +3. Write `ascending/r10m/zarr.json` with `spatial:shape` and `spatial:transform` +4. Create arrays: `vv/zarr.json`, `vh/zarr.json`, `border_mask/zarr.json` with `dimension_names: ["time", "y", "x"]` +5. Create coordinate arrays: `time/zarr.json` with `dimension_names: ["time"]`, etc. +5b. Create 1D spatial coordinate arrays (`x`, `y`) at every resolution level +6. Write first data shard: `vv/c/0/0/0` +7. Generate overviews for this timestep → write to r20m, r60m, ..., r720m + +**Append** (subsequent acquisitions): +1. Read current time dimension size +2. Resize time axis (increment shape[0]) +3. Write new shard: `vv/c/{new_index}/0/0` +4. Append coordinate values +5. Generate overviews for new timestep at all levels + +**Consolidation** (after all ingestions for a batch): +1. `zarr.consolidate_metadata(store_path, path=orbit_direction, zarr_format=3)` — orbit level +2. `zarr.consolidate_metadata(store_path, zarr_format=3)` — root level + +Note: Do NOT consolidate between individual ingestions — consolidated metadata caches array shapes and breaks subsequent resize+write operations. + +**CRITICAL — spatial:transform ordering**: S1Tiling GeoTIFFs use GDAL GeoTransform ordering `[c, a, b, f, d, e]`. The mini-spec requires Rasterio/Affine ordering `[a, b, c, d, e, f]`. The converter MUST apply the mapping: `spatial_transform = [GT(1), GT(2), GT(0), GT(4), GT(5), GT(3)]`. + +### 3.3 CLI Extension + +```bash +# Ingest a single acquisition +eopf-geozarr ingest-s1 \ + --vv /path/to/s1a_32TQM_vv_ASC_037_..._GammaNaughtRTC.tif \ + --vh /path/to/s1a_32TQM_vh_ASC_037_..._GammaNaughtRTC.tif \ + --mask /path/to/..._BorderMask.tif \ + --store s3://eopf-explorer/s1-grd-rtc-32TQM.zarr \ + --tile 32TQM \ + --orbit-dir ascending + +# Ingest conditions (once per tile/orbit) +eopf-geozarr ingest-s1-conditions \ + --lia /path/to/sin_LIA_32TQM_037.tif \ + --ia /path/to/IA_32TQM_037.tif \ + --gamma-area /path/to/GAMMA_AREA_32TQM_037.tif \ + --store s3://eopf-explorer/s1-grd-rtc-32TQM.zarr \ + --tile 32TQM \ + --orbit-dir ascending + +# Validate S1 store against mini-spec +eopf-geozarr validate-s1 s3://eopf-explorer/s1-grd-rtc-32TQM.zarr +``` + +--- + +## 4. Overview Generation + +Reuse existing downsampling logic with the same variable factor chain as S2: +- r10m → r20m (2×) → r60m (3×) → r120m (2×) → r360m (3×) → r720m (2×) + +Per-timestep generation: when a new acquisition is appended to r10m, generate overviews for that single timestep at all levels. Overviews are spatial only — no temporal resampling. + +Resampling methods: +- `vv`, `vh`: `"average"` (default in multiscales metadata) +- `border_mask`: `"nearest"` (binary mask, must not interpolate) + +--- + +## 5. STAC Registration + +New collection: `sentinel-1-grd-rtc-geozarr` + +Extensions: `sar`, `sat`, `zarr`, `render` + +Each store = one STAC item per tile. Temporal extent derived from sorted `time` coordinate. Updated on each ingest. + +--- + +## 6. Risks + +| Risk | Severity | Mitigation | Phase 0 Status | +|------|----------|------------|----------------| +| Zarr V3 append + sharding maturity | HIGH | Pin zarr-python, test heavily, serialize per tile | **MITIGATED** — resize+sharding verified with real 10980×10980 data, 3 timesteps appended successfully | +| Partial tile coverage per timestep | HIGH | border_mask in quality, coverage % in STAC | **MITIGATED** — real data shows 0.4%–99.8% coverage; border_mask ingestion validated | +| spatial:transform GDAL↔Rasterio ordering | MEDIUM | Explicit conversion in ingestion code, validation | **MITIGATED** — rasterio `src.transform` returns Affine ordering directly | +| Zarr V3 string dtype instability | MEDIUM | Use ` `UnstableSpecificationWarning: The data type (FixedLengthUTF32) does not have a Zarr V3 specification.` + +**Decision needed:** Use ` str: + """Classify an S1Tiling output file by its filename pattern.""" + name = path.name + if FNAME_PATTERN_MASK.match(name): + return "border_mask" + if FNAME_PATTERN_ORTHO.match(name): + if "_GammaNaughtRTC" in name: + return "gamma_naught_rtc" + if "_NormLim" in name: + return "sigma_normlim" + return "ortho" + if FNAME_PATTERN_LIA.match(name): + return "lia" + if FNAME_PATTERN_GAMMA_AREA.match(name): + return "gamma_area" + return "unknown" + + +def parse_filename_metadata(path: Path) -> dict: + """Extract metadata fields encoded in the filename.""" + name = path.name + result = {} + + m = FNAME_PATTERN_ORTHO.match(name) or FNAME_PATTERN_MASK.match(name) + if m: + result["flying_unit_code"] = f"s1{m.group(1)}" + result["tile_name"] = m.group(2) + result["polarisation"] = m.group(3).lower() + result["orbit_direction"] = m.group(4) + result["orbit"] = m.group(5) + result["acquisition_stamp"] = m.group(6) + return result + + m = FNAME_PATTERN_LIA.match(name) + if m: + result["lia_kind"] = m.group(1) + result["tile_name"] = m.group(2) + result["orbit"] = m.group(3) + return result + + m = FNAME_PATTERN_GAMMA_AREA.match(name) + if m: + result["tile_name"] = m.group(1) + result["orbit"] = m.group(2) + return result + + return result + + +def inspect_geotiff(path: Path) -> dict: + """Read all metadata from a GeoTIFF using rasterio.""" + with rasterio.open(str(path)) as src: + t = src.transform + info = { + "file": str(path), + "filename": path.name, + "file_type": classify_file(path), + "filename_metadata": parse_filename_metadata(path), + "rasterio_profile": { + "driver": src.driver, + "dtype": str(src.dtypes[0]) if src.dtypes else None, + "width": src.width, + "height": src.height, + "count": src.count, + "crs": str(src.crs) if src.crs else None, + "nodata": src.nodata, + }, + "transform": { + "affine": [t.a, t.b, t.c, t.d, t.e, t.f], + "pixel_size_x": abs(t.a), + "pixel_size_y": abs(t.e), + }, + "bounds": { + "left": src.bounds.left, + "bottom": src.bounds.bottom, + "right": src.bounds.right, + "top": src.bounds.top, + }, + "tags": dict(src.tags()), + "band_descriptions": list(src.descriptions) if src.descriptions else [], + "band_tags": {}, + } + # Per-band tags + for i in range(1, src.count + 1): + band_tags = src.tags(i) + if band_tags: + info["band_tags"][f"band_{i}"] = dict(band_tags) + + return info + + +def validate_tags(info: dict) -> list[str]: + """Validate that expected tags are present based on file type.""" + file_type = info["file_type"] + tags = info["tags"] + issues = [] + + if file_type == "unknown": + issues.append(f"Could not classify file: {info['filename']}") + return issues + + expected = { + "gamma_naught_rtc": EXPECTED_TAGS_ORTHO, + "sigma_normlim": EXPECTED_TAGS_ORTHO, + "ortho": EXPECTED_TAGS_ORTHO, + "border_mask": EXPECTED_TAGS_MASK, + "lia": EXPECTED_TAGS_LIA, + "gamma_area": EXPECTED_TAGS_GAMMA_AREA, + }.get(file_type, {}) + + # Check which expected tags are present/missing + for tag_name in expected: + if tag_name not in tags: + # Some tags are conditional — don't flag as errors + if tag_name in ("GAMMA_AREA_FILE", "ACQUISITION_DATETIME_2", "LIA_FILE"): + continue + issues.append(f"MISSING expected tag: {tag_name}") + + # Report unexpected tags (informational) + expected_names = set(expected.keys()) + for tag_name in tags: + if tag_name not in expected_names: + issues.append(f"EXTRA tag found: {tag_name} = {tags[tag_name]!r}") + + # Cross-validate filename metadata against tags + fname_meta = info["filename_metadata"] + if "flying_unit_code" in fname_meta and "FLYING_UNIT_CODE" in tags: + if fname_meta["flying_unit_code"] != tags["FLYING_UNIT_CODE"]: + issues.append( + f"MISMATCH flying_unit_code: filename={fname_meta['flying_unit_code']!r} " + f"vs tag={tags['FLYING_UNIT_CODE']!r}" + ) + if "orbit_direction" in fname_meta and "ORBIT_DIRECTION" in tags: + # Filename may use ASC/DES, tag may use ASCENDING/DESCENDING or vice versa + fname_dir = fname_meta["orbit_direction"].upper() + tag_dir = tags["ORBIT_DIRECTION"].upper() + if fname_dir not in tag_dir and tag_dir not in fname_dir: + issues.append(f"MISMATCH orbit_direction: filename={fname_dir!r} vs tag={tag_dir!r}") + + return issues + + +def print_report(info: dict, validate: bool = False) -> None: + """Print a human-readable report for one file.""" + print(f"\n{'=' * 80}") + print(f"FILE: {info['filename']}") + print(f"TYPE: {info['file_type']}") + print(f"PATH: {info['file']}") + print(f"{'=' * 80}") + + print("\n--- Rasterio Profile ---") + for k, v in info["rasterio_profile"].items(): + print(f" {k}: {v}") + + print("\n--- Transform ---") + print(f" Affine [a,b,c,d,e,f]: {info['transform']['affine']}") + print( + f" Pixel size: {info['transform']['pixel_size_x']}m x {info['transform']['pixel_size_y']}m" + ) + + print("\n--- Bounds ---") + for k, v in info["bounds"].items(): + print(f" {k}: {v}") + + print("\n--- Filename Metadata ---") + for k, v in info["filename_metadata"].items(): + print(f" {k}: {v}") + + print(f"\n--- GeoTIFF Tags ({len(info['tags'])} total) ---") + for k in sorted(info["tags"].keys()): + print(f" {k}: {info['tags'][k]!r}") + + if info["band_descriptions"]: + print("\n--- Band Descriptions ---") + for i, desc in enumerate(info["band_descriptions"], 1): + print(f" Band {i}: {desc}") + + if info["band_tags"]: + print("\n--- Band Tags ---") + for band, tags in info["band_tags"].items(): + print(f" {band}:") + for k, v in tags.items(): + print(f" {k}: {v!r}") + + if validate: + issues = validate_tags(info) + print(f"\n--- Validation ({len(issues)} issues) ---") + if not issues: + print(" ALL OK: all expected tags present, no mismatches") + else: + for issue in issues: + print(f" {'WARNING' if 'EXTRA' in issue else 'ERROR'}: {issue}") + + # Critical tags for our ingestion pipeline + print("\n--- Ingestion-Critical Tags ---") + critical = [ + "ACQUISITION_DATETIME", + "ORBIT_NUMBER", + "RELATIVE_ORBIT_NUMBER", + "FLYING_UNIT_CODE", + "ORBIT_DIRECTION", + "POLARIZATION", + "S2_TILE_CORRESPONDING_CODE", + "CALIBRATION", + "IMAGE_TYPE", + ] + for tag in critical: + value = info["tags"].get(tag) + status = "PRESENT" if value is not None else "MISSING" + print(f" [{status:>7}] {tag}: {value!r}") + + +def main(): + parser = argparse.ArgumentParser( + description="Inspect S1Tiling GeoTIFF outputs for metadata validation" + ) + parser.add_argument( + "path", + type=Path, + help="Path to a GeoTIFF file or directory containing GeoTIFFs", + ) + parser.add_argument( + "--validate", + action="store_true", + help="Validate tags against expected S1Tiling schema", + ) + parser.add_argument( + "--json", + action="store_true", + help="Output as JSON instead of human-readable format", + ) + args = parser.parse_args() + + if not args.path.exists(): + print(f"Error: path does not exist: {args.path}", file=sys.stderr) + sys.exit(1) + + # Collect files + if args.path.is_file(): + files = [args.path] + else: + files = sorted(args.path.rglob("*.tif")) + if not files: + print(f"No .tif files found in {args.path}", file=sys.stderr) + sys.exit(1) + + results = [] + for f in files: + info = inspect_geotiff(f) + results.append(info) + + if args.json: + print(json.dumps(results, indent=2, default=str)) + else: + print(f"Found {len(results)} GeoTIFF file(s)") + for info in results: + print_report(info, validate=args.validate) + + # Summary table + if len(results) > 1: + print(f"\n{'=' * 80}") + print("SUMMARY") + print(f"{'=' * 80}") + print(f" Files inspected: {len(results)}") + types = {} + for info in results: + ft = info["file_type"] + types[ft] = types.get(ft, 0) + 1 + for ft, count in sorted(types.items()): + print(f" {ft}: {count}") + + # Union of all tag keys seen + all_tags = set() + for info in results: + all_tags.update(info["tags"].keys()) + print(f"\n All tag keys across files ({len(all_tags)}):") + for tag in sorted(all_tags): + present_in = sum(1 for info in results if tag in info["tags"]) + print(f" {tag}: present in {present_in}/{len(results)} files") + + +if __name__ == "__main__": + main() diff --git a/analysis/s1_grd_rtc_prototype.py b/analysis/s1_grd_rtc_prototype.py new file mode 100644 index 00000000..10ebd671 --- /dev/null +++ b/analysis/s1_grd_rtc_prototype.py @@ -0,0 +1,632 @@ +""" +Phase 0 Prototype: S1Tiling GeoTIFF → GeoZarr V3 Store + +This script demonstrates the end-to-end GeoTIFF → Zarr conversion for S1 GRD RTC data. +It validates key technical assumptions: + +1. Zarr V3 array creation with sharding (zarr-python 3.1.1) +2. Time-axis resize/append model (shape[0] starts at 0, grows per acquisition) +3. GeoTIFF metadata extraction with rasterio (CRS, transform, custom tags) +4. GDAL → Rasterio/Affine spatial:transform conversion +5. GeoZarr Zarr Conventions (multiscales, proj:, spatial:) attribute structure +6. xarray compatibility for reading back the store +7. Overview generation with variable downsampling factors + +Run: python analysis/s1_grd_rtc_prototype.py +""" + +from __future__ import annotations + +import shutil +import tempfile +from pathlib import Path + +import numpy as np +import rasterio +import xarray as xr +import zarr +from rasterio.transform import from_bounds + +# ============================================================================= +# Constants: Zarr Conventions +# ============================================================================= + +MULTISCALES_UUID = "d35379db-88df-4056-af3a-620245f8e347" +GEO_PROJ_UUID = "f17cb550-5864-4468-aeb7-f3180cfb622f" +SPATIAL_UUID = "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4" + +ZARR_CONVENTIONS = [ + { + "uuid": MULTISCALES_UUID, + "schema_url": "https://raw.githubusercontent.com/zarr-conventions/multiscales/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-conventions/multiscales/blob/v1/README.md", + "name": "multiscales", + "description": "Multiscale layout of zarr datasets", + }, + { + "uuid": GEO_PROJ_UUID, + "schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md", + "name": "proj:", + "description": "Coordinate reference system information for geospatial data", + }, + { + "uuid": SPATIAL_UUID, + "schema_url": "https://raw.githubusercontent.com/zarr-conventions/spatial/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-conventions/spatial/blob/v1/README.md", + "name": "spatial:", + "description": "Spatial coordinate information", + }, +] + +# Overview chain: (level_name, parent_name, downsample_factor) +OVERVIEW_CHAIN = [ + ("r10m", None, 1), + ("r20m", "r10m", 2), + ("r60m", "r20m", 3), + ("r120m", "r60m", 2), + ("r360m", "r120m", 3), + ("r720m", "r360m", 2), +] + + +# ============================================================================= +# GeoTIFF Helpers +# ============================================================================= + + +def create_test_geotiff( + path: str | Path, + data: np.ndarray, + crs: str = "EPSG:32633", + transform: rasterio.transform.Affine | None = None, + tags: dict[str, str] | None = None, +) -> None: + """Write a single-band GeoTIFF with optional metadata tags.""" + with rasterio.open( + str(path), + "w", + driver="GTiff", + height=data.shape[0], + width=data.shape[1], + count=1, + dtype=data.dtype, + crs=crs, + transform=transform, + ) as dst: + if tags: + dst.update_tags(**tags) + dst.write(data, 1) + + +def extract_geotiff_metadata(path: str | Path) -> dict: + """Extract CRS, transform, bounds, and custom tags from a GeoTIFF.""" + with rasterio.open(str(path)) as src: + tags = src.tags() + t = src.transform + # spatial:transform in Rasterio/Affine ordering [a, b, c, d, e, f] + # This IS the native rasterio transform — no GDAL conversion needed + # when reading with rasterio (rasterio already returns Affine ordering). + spatial_transform = [t.a, t.b, t.c, t.d, t.e, t.f] + return { + "crs": str(src.crs), + "spatial_transform": spatial_transform, + "shape": list(src.shape), + "bounds": [src.bounds.left, src.bounds.bottom, src.bounds.right, src.bounds.top], + "datetime": tags.get("ACQUISITION_DATETIME"), + "absolute_orbit": int(tags.get("ORBIT_NUMBER", 0)), + "relative_orbit": int(tags.get("RELATIVE_ORBIT_NUMBER", 0)), + "platform": tags.get("FLYING_UNIT_CODE", ""), + } + + +# ============================================================================= +# Zarr Store Creation +# ============================================================================= + + +def compute_multiscales_layout( + native_shape: list[int], + native_transform: list[float], +) -> list[dict]: + """Build the multiscales layout array for all resolution levels.""" + layout: list[dict] = [] + current_shape = native_shape[:] + current_transform = native_transform[:] + + for level_name, parent_name, factor in OVERVIEW_CHAIN: + if parent_name is not None: + # Compute downsampled shape (ceiling division to match S2 convention) + current_shape = [ + int(np.ceil(current_shape[0] / factor)), + int(np.ceil(current_shape[1] / factor)), + ] + # Update transform: scale pixel size by factor, keep origin + current_transform = [ + current_transform[0] * factor, # a: pixel width + current_transform[1], # b: rotation (0) + current_transform[2], # c: x origin + current_transform[3], # d: rotation (0) + current_transform[4] * factor, # e: pixel height (negative) + current_transform[5], # f: y origin + ] + + entry: dict = { + "asset": level_name, + "spatial:shape": current_shape[:], + "spatial:transform": current_transform[:], + } + if parent_name is None: + entry["transform"] = {"scale": [1.0, 1.0]} + else: + entry["derived_from"] = parent_name + entry["transform"] = { + "scale": [float(factor), float(factor)], + "translation": [0.0, 0.0], + } + + layout.append(entry) + + return layout + + +def create_s1_store( + store_path: str | Path, + orbit_direction: str, + meta: dict, +) -> None: + """Create a new S1 GRD RTC Zarr V3 store with full conventions metadata.""" + height, width = meta["shape"] + + root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + orbit_group = root.create_group(orbit_direction) + + # Build full multiscales layout + 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"], + } + ) + + # Create each resolution level + for level_entry in layout: + level_name = level_entry["asset"] + level_shape = level_entry["spatial:shape"] + level_h, level_w = level_shape + + level_group = orbit_group.create_group(level_name) + level_group.attrs.update( + { + "spatial:shape": level_shape, + "spatial:transform": level_entry["spatial:transform"], + } + ) + + inner_chunks = (1, min(512, level_h), min(512, level_w)) + shard_shape = (1, level_h, level_w) + + # Data arrays + 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"], + ) + + # 1D spatial coordinate arrays (x and y) + # Compute from the level's spatial:transform [a, b, c, d, e, f] + lvl_transform = level_entry["spatial:transform"] + pixel_w = lvl_transform[0] # a: pixel width + x_origin = lvl_transform[2] # c: x origin (left edge) + pixel_h = lvl_transform[4] # e: pixel height (negative) + y_origin = lvl_transform[5] # f: y origin (top edge) + + x_coords = np.linspace( + x_origin, x_origin + level_w * pixel_w, level_w, endpoint=False, dtype="float64" + ) + y_coords = np.linspace( + y_origin, y_origin + level_h * pixel_h, level_h, endpoint=False, dtype="float64" + ) + + x_arr = level_group.create_array( + "x", + data=x_coords, + chunks=(level_w,), + fill_value=float("nan"), + dimension_names=["x"], + ) + x_arr.attrs.update( + { + "units": "m", + "long_name": "x coordinate of projection", + "standard_name": "projection_x_coordinate", + "_ARRAY_DIMENSIONS": ["x"], + } + ) + + y_arr = level_group.create_array( + "y", + data=y_coords, + chunks=(level_h,), + fill_value=float("nan"), + dimension_names=["y"], + ) + y_arr.attrs.update( + { + "units": "m", + "long_name": "y coordinate of projection", + "standard_name": "projection_y_coordinate", + "_ARRAY_DIMENSIONS": ["y"], + } + ) + + # Coordinate variables at native resolution only + r10m = orbit_group["r10m"] + for name, dtype, fill in [ + ("time", "int64", 0), + ("absolute_orbit", "int32", 0), + ("relative_orbit", "int32", 0), + ]: + r10m.create_array( + name, + shape=(0,), + dtype=dtype, + chunks=(512,), + fill_value=fill, + dimension_names=["time"], + ) + # Platform uses variable-length bytes to avoid Zarr V3 string dtype instability + r10m.create_array( + "platform", + shape=(0,), + dtype=" np.ndarray: + """Downsample a 2D array by the given factor.""" + h, w = data.shape + new_h = int(np.ceil(h / factor)) + new_w = int(np.ceil(w / factor)) + + if method == "nearest": + return data[::factor, ::factor][:new_h, :new_w] + + # Average: use block mean, handling edge blocks with padding + pad_h = new_h * factor - h + pad_w = new_w * factor - w + if pad_h > 0 or pad_w > 0: + padded = np.pad(data, ((0, pad_h), (0, pad_w)), mode="edge") + else: + padded = data + + reshaped = padded.reshape(new_h, factor, new_w, factor) + if np.issubdtype(data.dtype, np.floating): + return np.nanmean(reshaped, axis=(1, 3)).astype(data.dtype) + return reshaped.mean(axis=(1, 3)).astype(data.dtype) + + +# ============================================================================= +# Ingestion +# ============================================================================= + + +def ingest_acquisition( + store_path: str | Path, + orbit_direction: str, + vv_path: str | Path, + vh_path: str | Path, + mask_path: str | Path, + meta: dict, +) -> int: + """Append one acquisition to the store, including overviews.""" + root = zarr.open_group(str(store_path), mode="r+", zarr_format=3) + orbit = root[orbit_direction] + + # Read GeoTIFF 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(mask_path)) as src: + mask_data = src.read(1).astype(np.uint8) + + # Determine new time index + r10m = orbit["r10m"] + current_size = r10m["vv"].shape[0] + new_size = current_size + 1 + + # Write native resolution + all overview levels + data_by_level = {"r10m": (vv_data, vh_data, mask_data)} + + # Generate overviews + prev_vv, prev_vh, prev_mask = vv_data, vh_data, mask_data + for level_name, _, factor in OVERVIEW_CHAIN[1:]: + prev_vv = downsample_2d(prev_vv, factor, "average") + prev_vh = downsample_2d(prev_vh, factor, "average") + prev_mask = downsample_2d(prev_mask, factor, "nearest") + data_by_level[level_name] = (prev_vv, prev_vh, prev_mask) + + # Write to each level + for level_name, (vv_lev, vh_lev, mask_lev) in data_by_level.items(): + level = orbit[level_name] + h, w = vv_lev.shape + + level["vv"].resize((new_size, h, w)) + level["vh"].resize((new_size, h, w)) + level["border_mask"].resize((new_size, h, w)) + + level["vv"][current_size, :, :] = vv_lev + level["vh"][current_size, :, :] = vh_lev + level["border_mask"][current_size, :, :] = mask_lev + + # Append coordinate variables at native resolution + for coord_name in ["time", "absolute_orbit", "relative_orbit", "platform"]: + r10m[coord_name].resize((new_size,)) + + dt_ns = ( + np.datetime64(meta["datetime"].replace("Z", "")).astype("datetime64[ns]").astype(np.int64) + ) + r10m["time"][current_size] = dt_ns + r10m["absolute_orbit"][current_size] = meta["absolute_orbit"] + r10m["relative_orbit"][current_size] = meta["relative_orbit"] + r10m["platform"][current_size] = meta["platform"] + + return current_size + + +def consolidate_store( + store_path: str | Path, + orbit_direction: str, +) -> None: + """Consolidate metadata at orbit direction and root levels.""" + zarr.consolidate_metadata(str(store_path), path=orbit_direction, zarr_format=3) + zarr.consolidate_metadata(str(store_path), zarr_format=3) + + +# ============================================================================= +# Validation +# ============================================================================= + + +def validate_store(store_path: str | Path) -> None: + """Validate store structure and metadata.""" + root = zarr.open_group(str(store_path), mode="r", zarr_format=3) + errors: list[str] = [] + + # Check root-level consolidated metadata + root_meta = root.metadata + if root_meta.consolidated_metadata is None: + errors.append("root: missing consolidated_metadata") + + for orbit_dir in ["ascending", "descending"]: + if orbit_dir not in root: + continue + orbit = root[orbit_dir] + attrs = dict(orbit.attrs) + + # Check orbit-level consolidated metadata + orbit_meta = orbit.metadata + if orbit_meta.consolidated_metadata is None: + errors.append(f"{orbit_dir}: missing consolidated_metadata") + + # Check zarr_conventions + if "zarr_conventions" not in attrs: + errors.append(f"{orbit_dir}: missing zarr_conventions") + else: + conv_names = {c["name"] for c in attrs["zarr_conventions"]} + for required in ["multiscales", "proj:", "spatial:"]: + if required not in conv_names: + errors.append(f"{orbit_dir}: missing convention {required}") + + # Check proj:code + if "proj:code" not in attrs: + errors.append(f"{orbit_dir}: missing proj:code") + + # Check spatial:dimensions + if "spatial:dimensions" not in attrs: + errors.append(f"{orbit_dir}: missing spatial:dimensions") + + # Check multiscales layout + multiscales = attrs.get("multiscales", {}) + layout = multiscales.get("layout", []) + if not layout: + errors.append(f"{orbit_dir}: empty multiscales layout") + for entry in layout: + asset = entry.get("asset", "?") + if asset not in orbit: + errors.append(f"{orbit_dir}: layout references missing group {asset}") + if "spatial:transform" not in entry: + errors.append(f"{orbit_dir}/{asset}: missing spatial:transform in layout") + + # Check array dimension_names and coordinate arrays + for level_name in orbit.keys(): + level = orbit[level_name] + if not isinstance(level, zarr.Group): + continue + for arr_name in ["vv", "vh", "border_mask"]: + if arr_name in level: + arr = level[arr_name] + dim_names = arr.metadata.dimension_names + if dim_names != ("time", "y", "x"): + errors.append( + f"{orbit_dir}/{level_name}/{arr_name}: " + f"expected dimension_names ('time', 'y', 'x'), got {dim_names}" + ) + # Check x and y coordinate arrays exist + for coord_name in ["x", "y"]: + if coord_name not in level: + errors.append( + f"{orbit_dir}/{level_name}: missing {coord_name} coordinate array" + ) + else: + coord_arr = level[coord_name] + if len(coord_arr.shape) != 1: + errors.append( + f"{orbit_dir}/{level_name}/{coord_name}: " + f"expected 1D, got shape {coord_arr.shape}" + ) + + if errors: + print("VALIDATION ERRORS:") + for e in errors: + print(f" - {e}") + else: + print("VALIDATION PASSED") + + +# ============================================================================= +# Main: End-to-End Test +# ============================================================================= + + +def main() -> None: + tmpdir = Path(tempfile.mkdtemp()) + store_path = tmpdir / "s1-grd-rtc-32TQM.zarr" + + try: + SIZE = 256 + xmin, ymin, xmax, ymax = 500000.0, 4997440.0, 502560.0, 5000000.0 + transform = from_bounds(xmin, ymin, xmax, ymax, SIZE, SIZE) + + # --- Create synthetic acquisitions --- + acq1_tags = { + "ACQUISITION_DATETIME": "2023-01-15T06:12:34Z", + "ORBIT_NUMBER": "47001", + "RELATIVE_ORBIT_NUMBER": "037", + "FLYING_UNIT_CODE": "S1A", + } + acq2_tags = { + "ACQUISITION_DATETIME": "2023-01-27T06:12:35Z", + "ORBIT_NUMBER": "47177", + "RELATIVE_ORBIT_NUMBER": "037", + "FLYING_UNIT_CODE": "S1A", + } + + np.random.seed(42) + for suffix, tags in [("acq1", acq1_tags), ("acq2", acq2_tags)]: + for pol in ["vv", "vh"]: + create_test_geotiff( + tmpdir / f"{pol}_{suffix}.tif", + np.random.uniform(0, 1, (SIZE, SIZE)).astype(np.float32), + transform=transform, + tags=tags, + ) + create_test_geotiff( + tmpdir / f"mask_{suffix}.tif", + np.ones((SIZE, SIZE), dtype=np.float32), + transform=transform, + tags=tags, + ) + + # --- Create store --- + meta1 = extract_geotiff_metadata(tmpdir / "vv_acq1.tif") + print(f"[1/6] Creating store: CRS={meta1['crs']}, shape={meta1['shape']}") + create_s1_store(store_path, "ascending", meta1) + print(" Store created with 6 resolution levels") + + # --- Ingest acquisition 1 --- + idx = ingest_acquisition( + store_path, + "ascending", + tmpdir / "vv_acq1.tif", + tmpdir / "vh_acq1.tif", + tmpdir / "mask_acq1.tif", + meta1, + ) + print(f"[2/6] Ingested acquisition 1 at time_index={idx}") + + # --- Ingest acquisition 2 (append) --- + meta2 = extract_geotiff_metadata(tmpdir / "vv_acq2.tif") + idx2 = ingest_acquisition( + store_path, + "ascending", + tmpdir / "vv_acq2.tif", + tmpdir / "vh_acq2.tif", + tmpdir / "mask_acq2.tif", + meta2, + ) + print(f"[3/6] Ingested acquisition 2 at time_index={idx2}") + + # --- Consolidate metadata --- + consolidate_store(store_path, "ascending") + print("[3.5/6] Consolidated metadata at orbit and root levels") + + # --- Validate store --- + print("[4/6] Validating store structure...") + validate_store(store_path) + + # --- Verify data integrity --- + root = zarr.open_group(str(store_path), mode="r", zarr_format=3) + r10m = root["ascending/r10m"] + assert r10m["vv"].shape == (2, SIZE, SIZE), f"Unexpected VV shape: {r10m['vv'].shape}" + assert r10m["time"].shape == (2,), f"Unexpected time shape: {r10m['time'].shape}" + print("[5/6] Data integrity verified (2 timesteps, all arrays consistent)") + + # --- Verify overview levels --- + expected_shapes = { + "r10m": (2, 256, 256), + "r20m": (2, 128, 128), + "r60m": (2, 43, 43), + "r120m": (2, 22, 22), + "r360m": (2, 8, 8), + "r720m": (2, 4, 4), + } + for level, expected in expected_shapes.items(): + actual = root[f"ascending/{level}/vv"].shape + assert actual == expected, f"{level}: expected {expected}, got {actual}" + print("[6/6] Overview levels verified:") + for level, shape in expected_shapes.items(): + print(f" {level}: {shape[1]}×{shape[2]} ({shape[0]} timesteps)") + + # --- Read with xarray --- + print() + ds = xr.open_zarr(str(store_path / "ascending" / "r10m"), zarr_format=3, consolidated=False) + print(f"xarray read OK: variables={list(ds.data_vars)}, dims={dict(ds.sizes)}") + + print() + print("=" * 60) + print("PHASE 0 PROTOTYPE: ALL CHECKS PASSED") + print("=" * 60) + print() + print("Validated:") + print(" ✓ Zarr V3 + sharding (zarr-python 3.1.1)") + print(" ✓ Time-axis resize/append model") + print(" ✓ GeoTIFF metadata extraction (rasterio)") + print(" ✓ GeoZarr conventions (multiscales, proj:, spatial:)") + print(" ✓ 6-level overview pyramid (2x/3x chain)") + print(" ✓ xarray interoperability") + print(" ✓ Data integrity across append operations") + + finally: + shutil.rmtree(tmpdir) + + +if __name__ == "__main__": + main() diff --git a/analysis/s1_real_geotiff_to_zarr.py b/analysis/s1_real_geotiff_to_zarr.py new file mode 100644 index 00000000..dd7328dc --- /dev/null +++ b/analysis/s1_real_geotiff_to_zarr.py @@ -0,0 +1,702 @@ +""" +Phase 0 — Real GeoTIFF → GeoZarr V3 Conversion Test + +Reads actual S1Tiling γ0T RTC GeoTIFFs produced by the Docker pipeline +and converts them into a Zarr V3 store following the implementation plan. + +Three acquisitions (orbits 008, 037, 110) × two polarisations (VV, VH) ++ border masks. Full 10980×10980 at 10m resolution, EPSG:32631. + +Usage: + python analysis/s1_real_geotiff_to_zarr.py \ + --input-dir ~/Downloads/s1tiling_test/data_out/31TCH \ + --gamma-area-dir ~/Downloads/s1tiling_test/data_gamma_area \ + --output-dir ~/Downloads/s1tiling_test/zarr_test +""" + +from __future__ import annotations + +import argparse +import re +import sys +import time +from pathlib import Path + +import numpy as np +import rasterio +import xarray as xr +import zarr + +# ============================================================================= +# Constants: Zarr Conventions (same as prototype) +# ============================================================================= + +MULTISCALES_UUID = "d35379db-88df-4056-af3a-620245f8e347" +GEO_PROJ_UUID = "f17cb550-5864-4468-aeb7-f3180cfb622f" +SPATIAL_UUID = "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4" + +ZARR_CONVENTIONS = [ + { + "uuid": MULTISCALES_UUID, + "schema_url": "https://raw.githubusercontent.com/zarr-conventions/multiscales/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-conventions/multiscales/blob/v1/README.md", + "name": "multiscales", + "description": "Multiscale layout of zarr datasets", + }, + { + "uuid": GEO_PROJ_UUID, + "schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md", + "name": "proj:", + "description": "Coordinate reference system information for geospatial data", + }, + { + "uuid": SPATIAL_UUID, + "schema_url": "https://raw.githubusercontent.com/zarr-conventions/spatial/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-conventions/spatial/blob/v1/README.md", + "name": "spatial:", + "description": "Spatial coordinate information", + }, +] + +# Overview chain: (level_name, parent_name, downsample_factor) +OVERVIEW_CHAIN = [ + ("r10m", None, 1), + ("r20m", "r10m", 2), + ("r60m", "r20m", 3), + ("r120m", "r60m", 2), + ("r360m", "r120m", 3), + ("r720m", "r360m", 2), +] + + +def best_chunk_size(dim: int, max_chunk: int = 512) -> int: + """Find the largest divisor of *dim* that is <= *max_chunk*. + + Zarr sharding requires inner chunks to divide the shard evenly. + """ + if dim <= max_chunk: + return dim + for d in range(max_chunk, 0, -1): + if dim % d == 0: + return d + return 1 + + +# Filename pattern for S1Tiling outputs +# e.g. s1a_31TCH_vv_DES_110_20250205t060110_GammaNaughtRTC.tif +S1TILING_PATTERN = re.compile( + r"(?Ps1[abc])_" + r"(?P[0-9]{2}[A-Z]{3})_" + r"(?Pvv|vh)_" + r"(?PASC|DES)_" + r"(?P\d{3})_" + r"(?P\d{8}t\d{6})_" + r"(?PGammaNaughtRTC)" + r"(?P_BorderMask)?\.tif$" +) + + +# ============================================================================= +# GeoTIFF metadata extraction +# ============================================================================= + + +def extract_geotiff_metadata(path: Path) -> dict: + """Extract CRS, transform, bounds, and custom tags from a real S1Tiling GeoTIFF.""" + with 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] + return { + "crs": str(src.crs), + "spatial_transform": spatial_transform, + "shape": [src.height, src.width], + "bounds": [src.bounds.left, src.bounds.bottom, src.bounds.right, src.bounds.top], + "datetime": tags.get("ACQUISITION_DATETIME", ""), + "absolute_orbit": int(tags.get("ORBIT_NUMBER", "0")), + "relative_orbit": int(tags.get("RELATIVE_ORBIT_NUMBER", "0")), + "platform": tags.get("FLYING_UNIT_CODE", ""), + "calibration": tags.get("CALIBRATION", ""), + "input_s1_images": tags.get("INPUT_S1_IMAGES", ""), + "tags": tags, + } + + +# ============================================================================= +# Group S1Tiling output files into acquisitions +# ============================================================================= + + +def discover_acquisitions(input_dir: Path) -> list[dict]: + """Parse filenames and group into acquisition bundles (vv, vh, masks).""" + files = sorted(input_dir.glob("*.tif")) + # Group by (platform, tile, orbit_dir, rel_orbit, acq_stamp) + groups: dict[tuple, dict] = {} + + for f in files: + m = S1TILING_PATTERN.match(f.name) + if not m: + print(f" SKIP: {f.name} (doesn't match pattern)") + continue + + key = ( + m.group("platform"), + m.group("tile"), + m.group("orbit_dir"), + m.group("rel_orbit"), + m.group("acq_stamp"), + ) + + if key not in groups: + groups[key] = { + "platform": m.group("platform"), + "tile": m.group("tile"), + "orbit_dir": m.group("orbit_dir"), + "rel_orbit": m.group("rel_orbit"), + "acq_stamp": m.group("acq_stamp"), + } + + pol = m.group("pol") + is_mask = m.group("mask") is not None + + if is_mask: + groups[key][f"{pol}_mask"] = f + else: + groups[key][pol] = f + + # Validate completeness + acquisitions = [] + for key, acq in sorted(groups.items()): + missing = [k for k in ("vv", "vh", "vv_mask", "vh_mask") if k not in acq] + if missing: + print(f" WARNING: Acquisition {key} missing: {missing}") + acquisitions.append(acq) + + return acquisitions + + +# ============================================================================= +# Multiscales layout +# ============================================================================= + + +def compute_multiscales_layout( + native_shape: list[int], + native_transform: list[float], +) -> list[dict]: + """Build the multiscales layout array for all resolution levels.""" + layout: list[dict] = [] + current_shape = native_shape[:] + current_transform = native_transform[:] + + for level_name, parent_name, factor in OVERVIEW_CHAIN: + if parent_name is not None: + current_shape = [ + int(np.ceil(current_shape[0] / factor)), + int(np.ceil(current_shape[1] / factor)), + ] + current_transform = [ + current_transform[0] * factor, + current_transform[1], + current_transform[2], + current_transform[3], + current_transform[4] * factor, + current_transform[5], + ] + + entry: dict = { + "asset": level_name, + "spatial:shape": current_shape[:], + "spatial:transform": current_transform[:], + } + if parent_name is None: + entry["transform"] = {"scale": [1.0, 1.0]} + else: + entry["derived_from"] = parent_name + entry["transform"] = { + "scale": [float(factor), float(factor)], + "translation": [0.0, 0.0], + } + layout.append(entry) + + return layout + + +# ============================================================================= +# Store creation +# ============================================================================= + + +def create_s1_store( + store_path: Path, + orbit_direction: str, + meta: dict, +) -> None: + """Create a new S1 GRD RTC Zarr V3 store with full conventions metadata.""" + height, width = meta["shape"] + + root = zarr.open_group(str(store_path), mode="w", zarr_format=3) + 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"], + } + ) + + # Create each resolution level + 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, best_chunk_size(level_h), best_chunk_size(level_w)) + 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"], + ) + + # Coordinate variables at native resolution only + r10m = orbit_group["r10m"] + for name, dtype, fill in [ + ("time", "int64", 0), + ("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=" np.ndarray: + """Downsample a 2D array by the given factor.""" + h, w = data.shape + new_h = int(np.ceil(h / factor)) + new_w = int(np.ceil(w / factor)) + + if method == "nearest": + return data[::factor, ::factor][:new_h, :new_w] + + # Average with edge padding + pad_h = new_h * factor - h + pad_w = new_w * factor - w + if pad_h > 0 or pad_w > 0: + padded = np.pad(data, ((0, pad_h), (0, pad_w)), mode="edge") + else: + padded = data + + reshaped = padded.reshape(new_h, factor, new_w, factor) + if np.issubdtype(data.dtype, np.floating): + return np.nanmean(reshaped, axis=(1, 3)).astype(data.dtype) + return reshaped.mean(axis=(1, 3)).astype(data.dtype) + + +# ============================================================================= +# Ingestion: append one acquisition +# ============================================================================= + + +def ingest_acquisition( + store_path: Path, + orbit_direction: str, + vv_path: Path, + vh_path: Path, + vv_mask_path: Path, + vh_mask_path: Path, + meta: dict, +) -> int: + """Append one acquisition to the store, including overviews. + + Uses the VV border mask (S1Tiling produces per-pol masks; they differ + slightly because each pol band has different no-data fringe patterns). + For our pipeline, we use the VV mask as the primary border mask. + """ + root = zarr.open_group(str(store_path), mode="r+", zarr_format=3) + orbit = root[orbit_direction] + + # Read GeoTIFF data + t0 = time.time() + 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(vv_mask_path)) as src: + mask_data = src.read(1).astype(np.uint8) + read_time = time.time() - t0 + + print( + f" Read GeoTIFFs: {read_time:.1f}s " + f"(vv: {vv_data.dtype} min={np.nanmin(vv_data):.4f} max={np.nanmax(vv_data):.4f}, " + f"mask: unique={np.unique(mask_data).tolist()})" + ) + + # Determine new time index + r10m = orbit["r10m"] + current_size = r10m["vv"].shape[0] + new_size = current_size + 1 + + # Build data at all levels + data_by_level = {"r10m": (vv_data, vh_data, mask_data)} + + t0 = time.time() + prev_vv, prev_vh, prev_mask = vv_data, vh_data, mask_data + for level_name, _, factor in OVERVIEW_CHAIN[1:]: + prev_vv = downsample_2d(prev_vv, factor, "average") + prev_vh = downsample_2d(prev_vh, factor, "average") + prev_mask = downsample_2d(prev_mask, factor, "nearest") + data_by_level[level_name] = (prev_vv, prev_vh, prev_mask) + overview_time = time.time() - t0 + print(f" Overviews generated: {overview_time:.1f}s") + + # Write to each level + t0 = time.time() + for level_name, (vv_lev, vh_lev, mask_lev) in data_by_level.items(): + level = orbit[level_name] + h, w = vv_lev.shape + + level["vv"].resize((new_size, h, w)) + level["vh"].resize((new_size, h, w)) + level["border_mask"].resize((new_size, h, w)) + + level["vv"][current_size, :, :] = vv_lev + level["vh"][current_size, :, :] = vh_lev + level["border_mask"][current_size, :, :] = mask_lev + write_time = time.time() - t0 + print(f" Zarr write (all levels): {write_time:.1f}s") + + # Append coordinate variables + for coord_name in ["time", "absolute_orbit", "relative_orbit", "platform"]: + r10m[coord_name].resize((new_size,)) + + # Parse datetime — S1Tiling uses format: "2025:02:10T06:09:20Z" + dt_str = meta["datetime"] + # Normalise separators: "2025:02:10T06:09:20Z" → "2025-02-10T06:09:20" + dt_normalised = dt_str.replace("Z", "") + # Handle "2025:02:10" date format from S1Tiling + parts = dt_normalised.split("T") + if len(parts) == 2: + date_part = parts[0].replace(":", "-") + dt_normalised = f"{date_part}T{parts[1]}" + + dt_ns = np.datetime64(dt_normalised).astype("datetime64[ns]").astype(np.int64) + r10m["time"][current_size] = dt_ns + r10m["absolute_orbit"][current_size] = meta["absolute_orbit"] + r10m["relative_orbit"][current_size] = meta["relative_orbit"] + r10m["platform"][current_size] = meta["platform"] + + return current_size + + +# ============================================================================= +# Ingestion: gamma_area conditions +# ============================================================================= + + +def ingest_gamma_area( + store_path: Path, + orbit_direction: str, + gamma_area_files: list[Path], + meta: dict, +) -> None: + """Write gamma_area condition arrays (time-invariant, per orbit).""" + root = zarr.open_group(str(store_path), mode="r+", zarr_format=3) + orbit = root[orbit_direction] + + # Create conditions group if it doesn't exist + if "conditions" not in orbit: + conditions = orbit.create_group("conditions") + conditions.attrs.update( + { + "proj:code": meta["crs"], + "spatial:dimensions": ["Y", "X"], + "spatial:transform": meta["spatial_transform"], + } + ) + else: + conditions = orbit["conditions"] + + for ga_path in gamma_area_files: + # Extract orbit number from filename: GAMMA_AREA_31TCH_008.tif + m = re.search(r"GAMMA_AREA_\w+_(\d{3})\.tif$", ga_path.name) + if not m: + print(f" SKIP gamma area: {ga_path.name}") + continue + + orbit_num = m.group(1) + array_name = f"gamma_area_{orbit_num}" + + with rasterio.open(str(ga_path)) as src: + data = src.read(1) + + print( + f" Writing {array_name}: shape={data.shape}, " + f"dtype={data.dtype}, min={np.nanmin(data):.4f}, max={np.nanmax(data):.4f}" + ) + + h, w = data.shape + if array_name in conditions: + conditions[array_name][:, :] = data + else: + arr = conditions.create_array( + array_name, + shape=(h, w), + dtype="float32", + chunks=(min(512, h), min(512, w)), + compressors=zarr.codecs.BloscCodec(cname="zstd", clevel=5), + fill_value=float("nan"), + dimension_names=["Y", "X"], + ) + arr[:, :] = data + + +# ============================================================================= +# Validation +# ============================================================================= + + +def validate_and_report(store_path: Path, orbit_direction: str) -> None: + """Validate the store and print a comprehensive report.""" + print("\n" + "=" * 72) + print("VALIDATION REPORT") + print("=" * 72) + + root = zarr.open_group(str(store_path), mode="r", zarr_format=3) + orbit = root[orbit_direction] + attrs = dict(orbit.attrs) + + # 1. Zarr conventions + convs = attrs.get("zarr_conventions", []) + conv_names = {c["name"] for c in convs} + for required in ["multiscales", "proj:", "spatial:"]: + status = "OK" if required in conv_names else "MISSING" + print(f" Convention '{required}': {status}") + + # 2. proj:code + print(f" proj:code: {attrs.get('proj:code', 'MISSING')}") + print(f" spatial:dimensions: {attrs.get('spatial:dimensions', 'MISSING')}") + print(f" spatial:bbox: {attrs.get('spatial:bbox', 'MISSING')}") + + # 3. Multiscales layout + layout = attrs.get("multiscales", {}).get("layout", []) + print(f"\n Multiscales layout ({len(layout)} levels):") + for entry in layout: + name = entry["asset"] + shape = entry["spatial:shape"] + res = entry["spatial:transform"][0] + derived = entry.get("derived_from", "—") + scale = entry["transform"].get("scale", [1, 1]) + print(f" {name}: {shape[0]}×{shape[1]} @ {res}m (from {derived}, scale {scale})") + + # 4. Array shapes and time dimension + print("\n Data arrays at r10m:") + r10m = orbit["r10m"] + for arr_name in ["vv", "vh", "border_mask"]: + arr = r10m[arr_name] + print( + f" {arr_name}: shape={arr.shape}, dtype={arr.dtype}, " + f"dim_names={arr.metadata.dimension_names}" + ) + + # 5. Coordinate variables + print("\n Coordinate variables:") + for coord in ["time", "absolute_orbit", "relative_orbit", "platform"]: + arr = r10m[coord] + vals = arr[:] + print(f" {coord}: shape={arr.shape}, dtype={arr.dtype}, values={vals}") + + # Decode time values + time_arr = r10m["time"][:] + if len(time_arr) > 0: + datetimes = time_arr.astype("datetime64[ns]") + print(f" time (decoded): {[str(dt) for dt in datetimes]}") + + # 6. Overview consistency + print("\n Overview shapes:") + for level_name, _, _ in OVERVIEW_CHAIN: + if level_name in orbit: + vv = orbit[level_name]["vv"] + print(f" {level_name}: vv.shape={vv.shape}") + + # 7. Conditions + if "conditions" in orbit: + conditions = orbit["conditions"] + print("\n Conditions:") + cond_attrs = dict(conditions.attrs) + print(f" proj:code: {cond_attrs.get('proj:code', 'MISSING')}") + for name in conditions.keys(): + item = conditions[name] + if hasattr(item, "shape"): + print(f" {name}: shape={item.shape}, dtype={item.dtype}") + + # 8. xarray readback + print("\n xarray readback:") + try: + ds = xr.open_zarr( + str(store_path / orbit_direction / "r10m"), zarr_format=3, consolidated=False + ) + print(f" Dataset: {dict(ds.dims)}") + print(f" Variables: {list(ds.data_vars)}") + print(f" Coords: {list(ds.coords)}") + # Read a small sample to verify data integrity + sample = ds["vv"].isel(time=0, Y=slice(0, 3), X=slice(0, 3)).values + print(f" vv[0, :3, :3] = {sample}") + ds.close() + except Exception as e: + print(f" ERROR: {e}") + + # 9. Store size on disk + store_size = sum(f.stat().st_size for f in store_path.rglob("*") if f.is_file()) + print(f"\n Total store size: {store_size / 1e9:.2f} GB") + + print("\n" + "=" * 72) + + +# ============================================================================= +# Main +# ============================================================================= + + +def main() -> None: + parser = argparse.ArgumentParser( + description="Convert real S1Tiling GeoTIFFs to GeoZarr V3 store" + ) + parser.add_argument( + "--input-dir", required=True, help="Path to S1Tiling output directory (e.g. data_out/31TCH)" + ) + parser.add_argument("--gamma-area-dir", default=None, help="Path to gamma area maps directory") + parser.add_argument( + "--output-dir", required=True, help="Directory where Zarr store will be created" + ) + args = parser.parse_args() + + input_dir = Path(args.input_dir).expanduser() + output_dir = Path(args.output_dir).expanduser() + output_dir.mkdir(parents=True, exist_ok=True) + + print("=" * 72) + print("S1 GRD RTC — Real GeoTIFF → GeoZarr V3 Conversion") + print("=" * 72) + + # 1. Discover acquisitions + print(f"\nDiscovering acquisitions in {input_dir}...") + acquisitions = discover_acquisitions(input_dir) + print(f" Found {len(acquisitions)} acquisitions:") + for acq in acquisitions: + print( + f" {acq['platform']} orbit {acq['rel_orbit']} @ {acq['acq_stamp']} ({acq['orbit_dir']})" + ) + + if not acquisitions: + print("ERROR: No acquisitions found!") + sys.exit(1) + + # 2. Extract metadata from first VV file to initialise the store + first_acq = acquisitions[0] + meta = extract_geotiff_metadata(first_acq["vv"]) + print("\n Reference metadata:") + print(f" CRS: {meta['crs']}") + print(f" Shape: {meta['shape']}") + print(f" Transform: {meta['spatial_transform']}") + print(f" Bounds: {meta['bounds']}") + print(f" Calibration: {meta['calibration']}") + + # Determine orbit direction + orbit_dir_code = first_acq["orbit_dir"] + orbit_direction = "descending" if orbit_dir_code == "DES" else "ascending" + + # Tile ID from first acquisition + tile_id = first_acq["tile"] + store_path = output_dir / f"s1-grd-rtc-{tile_id}.zarr" + + # 3. Create the store + print(f"\nCreating store: {store_path}") + total_t0 = time.time() + create_s1_store(store_path, orbit_direction, meta) + print(" Store structure created.") + + # 4. Ingest each acquisition + for i, acq in enumerate(acquisitions): + print( + f"\n--- Ingesting acquisition {i + 1}/{len(acquisitions)}: " + f"{acq['platform']} orbit {acq['rel_orbit']} @ {acq['acq_stamp']} ---" + ) + + acq_meta = extract_geotiff_metadata(acq["vv"]) + t0 = time.time() + idx = ingest_acquisition( + store_path, + orbit_direction, + vv_path=acq["vv"], + vh_path=acq["vh"], + vv_mask_path=acq["vv_mask"], + vh_mask_path=acq["vh_mask"], + meta=acq_meta, + ) + acq_time = time.time() - t0 + print(f" → Time index {idx}, acquisition total: {acq_time:.1f}s") + + # 5. Ingest gamma_area conditions + if args.gamma_area_dir: + gamma_dir = Path(args.gamma_area_dir).expanduser() + gamma_files = sorted(gamma_dir.glob("GAMMA_AREA_*.tif")) + if gamma_files: + print(f"\n--- Ingesting {len(gamma_files)} gamma area condition(s) ---") + ingest_gamma_area(store_path, orbit_direction, gamma_files, meta) + + total_time = time.time() - total_t0 + print(f"\n Total conversion time: {total_time:.1f}s") + + # 6. Validate + validate_and_report(store_path, orbit_direction) + + +if __name__ == "__main__": + main() diff --git a/analysis/s1tiling_docker_instructions.md b/analysis/s1tiling_docker_instructions.md new file mode 100644 index 00000000..65108514 --- /dev/null +++ b/analysis/s1tiling_docker_instructions.md @@ -0,0 +1,406 @@ +# S1Tiling γ0T RTC Processing – Docker Instructions + +Run S1Tiling 1.4.0 via Docker to produce γ0T RTC-calibrated GeoTIFFs from a +Sentinel-1 GRD product downloaded from CDSE, orthorectified onto an S2 MGRS tile. + +> **Goal:** obtain *real* S1Tiling outputs to validate metadata/format +> assumptions for the GeoZarr ingestion pipeline. + +--- + +## Prerequisites + +| Requirement | Notes | +|---|---| +| Docker | `docker --version` ≥ 20 | +| Disk space | ~20 GB (S1 GRD ≈ 1.7 GB, DEM tiles ≈ 300 MB, outputs + tmp ≈ 15 GB) | +| RAM | ≥ 16 GB recommended (γ area estimation is RAM-greedy) | +| CDSE account | Register at | +| Internet | For S1 product download and DEM fetching | + +--- + +## 1 – Create working directory structure + +```bash +export S1T_WORKDIR=$HOME/Downloads/s1tiling_test +mkdir -p "$S1T_WORKDIR"/{data_out,data_raw,data_gamma_area,tmp,eof,config} +mkdir -p "$S1T_WORKDIR"/DEM/SRTM_30_hgt +``` + +--- + +## 2 – Configure EODAG for CDSE + +Create `$S1T_WORKDIR/config/eodag.yml`: + +```yaml +cop_dataspace: + priority: 1 + auth: + credentials: + username: "YOUR_CDSE_EMAIL" + password: "YOUR_CDSE_PASSWORD" +``` + +Replace the credentials with your CDSE (Copernicus Data Space Ecosystem) +account. EODAG will use `cop_dataspace` as the preferred provider to search +and download Sentinel-1 GRD products. The patch script fixes all the EODAG +4.0.0 incompatibilities so no other provider configuration is needed. + +> **Tip:** You can test your credentials at +> before running S1Tiling. + +--- + +## 3 – (Optional) Pre-download SRTM DEM tiles + +S1Tiling needs SRTM 30 m DEM tiles in `$S1T_WORKDIR/DEM/SRTM_30_hgt/`. + +If you already have them, copy/symlink. Otherwise S1Tiling can sometimes +auto-download, but it's more reliable to pre-stage them. + +For MGRS tile **31TCH** (south of France), 20 SRTM tiles are needed because +the S1 GRD swaths extend well beyond the MGRS tile footprint — the Gamma Area +computation (for RTC) needs DEM covering the full S1 acquisition geometry: + +```bash +cd "$S1T_WORKDIR/DEM/SRTM_30_hgt" +for tile in N41E002 N41E003 \ + N42E000 N42E001 N42E002 N42E003 N42W001 N42W002 N42W003 \ + N43E000 N43E001 N43E002 N43E003 N43E004 N43E005 N43W001 N43W002 N43W003 \ + N44W001 N44W002; do + lat="${tile:0:3}" + curl -sS -o "${tile}.hgt.gz" \ + "https://s3.amazonaws.com/elevation-tiles-prod/skadi/${lat}/${tile}.hgt.gz" + gunzip -f "${tile}.hgt.gz" +done +# Verify: should be 20 .hgt files (~25 MB each, ~500 MB total) +ls -1 *.hgt | wc -l +``` + +> **Why so many tiles?** The MGRS tile 31TCH only covers lat 42–43°N, +> lon 0–2°E, but the 3 overlapping S1 descending orbits (008, 037, 110) +> have swaths extending from ~41°N to ~44°N and from ~3°W to ~5°E. +> The AgglomerateDEM step in S1Tiling needs DEM for the full swath. + +> **Alternative tile:** Use any MGRS tile you prefer. Adjust DEM tiles and +> the config accordingly. Just pick a tile with Sentinel-1 coverage in the +> chosen date range. + +--- + +## 4 – Create the S1Tiling configuration file + +Create `$S1T_WORKDIR/config/S1GRD_RTC.cfg`: + +```ini +[Paths] +# Final γ0T RTC calibrated products +output : /data/data_out + +# Gamma Area maps directory +gamma_area : /data/data_gamma_area + +# Raw S1 products (downloaded here) +s1_images : /data/data_raw + +# Precise Orbit files (downloaded automatically) +eof_dir : /data/eof + +# Temporary files (can be large ~15 GB) +tmp : /tmp/s1tiling + +# DEM information +dem_dir : /MNT/SRTM_30_hgt +dem_info : SRTM 30m + +# Geoid (shipped inside the docker image — use the image's install path, NOT /data) +geoid_file : /opt/S1TilingEnv/lib/python3.10/site-packages/s1tiling/resources/Geoid/egm96.grd + +[DataSource] +# EODAG config (mounted in docker) +eodag_config : /eo_config/eodag.yml + +# Enable downloading from CDSE +download : True +nb_parallel_downloads : 2 + +# Region of interest: S2 MGRS tile(s) +# S1Tiling will download S1 GRD products overlapping this tile +roi_by_tiles : 31TCH + +# Platform filter (leave empty for both S1A and S1B) +platform_list : S1A + +# Polarisation +polarisation : VV VH + +# Orbit direction filter (optional — DES = descending) +orbit_direction : DES + +# Date range — keep it VERY SHORT (12 days = 1 revisit cycle) +# This minimises download volume. One date pair is enough for our test. +first_date : 2025-02-01 +last_date : 2025-02-14 + +[Processing] +# --- γ0T RTC calibration --- +calibration : gamma_naught_rtc + +# Noise removal +remove_thermal_noise : True +lower_signal_value : 1e-7 + +# Output resolution (meters) +output_spatial_resolution : 10. + +# Tiles to process +tiles : 31TCH + +# Orthorectification interpolation +orthorectification_interpolation_method : bco +orthorectification_gridspacing : 40 + +# DEM cache strategy +cache_dem_by : copy + +# γ area RTC specific +distribute_area : False +min_gamma_area : 1.0 +calibration_factor : 1.0 +output_nodata : False + +# Do NOT use resampled DEM (simpler, less RAM issue) +use_resampled_dem : False + +# Streaming: disable for gamma_area to avoid artefacts +disable_streaming.apply_gamma_area : True + +# Parallelism — for a single test, keep low +nb_parallel_processes : 1 +ram_per_process : 8192 +nb_otb_threads : 4 + +# Logging +mode : debug logging + +# Generate border masks (useful for our pipeline) +[Mask] +generate_border_mask : True + +[Quicklook] +generate : False + +[Metadata] +phase0_test : s1-grd-rtc-validation +producer : eopf-geozarr-phase0 +``` + +### Key options explained + +| Option | Value | Why | +|---|---|---| +| `calibration` | `gamma_naught_rtc` | This is the γ0T RTC pipeline | +| `roi_by_tiles` | `31TCH` | S2 MGRS tile (matches S1Tiling docs examples) | +| `first/last_date` | 12-day window | Minimises downloads — 1 revisit cycle | +| `platform_list` | `S1A` | Single platform to reduce data volume | +| `orbit_direction` | `DES` | Single direction to further reduce volume | +| `generate_border_mask` | `True` | Produces mask files we need for pipeline | + +--- + +## 5 – Pull the Docker image + +```bash +docker pull registry.orfeo-toolbox.org/s1-tiling/s1tiling:1.4.0-ubuntu-otb9.1.1 +``` + +Alternative registry: +```bash +docker pull cnes/s1tiling:1.4.0-ubuntu-otb9.1.1 +``` + +--- + +## 6 – Run S1Tiling (γ0T RTC) + +S1Tiling with `gamma_naught_rtc` calibration automatically computes the Gamma +Area maps and applies them. A single `S1Processor` invocation handles +everything. + +```bash +docker run --rm \ + -v "$S1T_WORKDIR"/DEM:/MNT \ + -v "$S1T_WORKDIR":/data \ + -v "$S1T_WORKDIR"/config:/eo_config \ + -v "$(dirname "$0")"/../analysis:/patch \ + --entrypoint bash \ + registry.orfeo-toolbox.org/s1-tiling/s1tiling:1.4.0-ubuntu-otb9.1.1 \ + -c 'python3 /patch/s1tiling_eodag4_patch.py && S1Processor /data/config/S1GRD_RTC.cfg' +``` + +If running from a different directory, replace the `-v` mount for `/patch` +with the absolute path to the `analysis/` folder: + +```bash +docker run --rm \ + -v "$S1T_WORKDIR"/DEM:/MNT \ + -v "$S1T_WORKDIR":/data \ + -v "$S1T_WORKDIR"/config:/eo_config \ + -v /home/emathot/Workspace/eopf-explorer/data-model/analysis:/patch \ + --entrypoint bash \ + registry.orfeo-toolbox.org/s1-tiling/s1tiling:1.4.0-ubuntu-otb9.1.1 \ + -c 'python3 /patch/s1tiling_eodag4_patch.py && S1Processor /data/config/S1GRD_RTC.cfg' +``` + +> **Bug workaround:** S1Tiling 1.4.0 ships EODAG 4.0.0 which has five +> breaking changes that prevent it from working with `cop_dataspace`: +> +> 1. `productType` kwarg to `dag.search()` was renamed to `collection`; +> having both causes cop_dataspace to fail silently → falls back to peps +> 2. Product properties now use STAC names (`sat:orbit_state`, `platform`, etc.) +> instead of legacy names (`orbitDirection`, `platformSerialIdentifier`, etc.) +> 3. `cop_dataspace` OData v4 API rejects `polarizationChannels` and `sensorMode` +> 4. `cop_dataspace` requires UPPERCASE orbit direction (`"DESCENDING"` not `"descending"`) +> 5. `relativeOrbitNumber` search param silently returns 0 results on cop_dataspace +> +> The patch script `s1tiling_eodag4_patch.py` fixes all five issues. +> The container is `--rm` so nothing persists. + +**What this does:** +1. Downloads S1A GRD products from CDSE (via EODAG) for the date range +2. Downloads precise orbit (EOF) files +3. Calibrates (σ0 internally), cuts, and orthorectifies onto the 31TCH MGRS grid +4. Computes the Gamma Area map for the orbit +5. Applies the γ0T RTC correction +6. Concatenates half-tiles into final products +7. Generates border masks + +**Expected runtime:** 30 min – 2 hours depending on network speed and host CPU. + +### Troubleshooting + +- **Download failures:** CDSE can be slow or rate-limited. Re-run the same + command — S1Tiling caches intermediary results. Already-completed steps are + skipped. +- **RAM issues:** If the gamma area computation fails with OOM, try setting + `use_resampled_dem : True` with `resample_dem_factor_x : 2.0` / + `resample_dem_factor_y : 2.0` in the config, or increase `ram_per_process`. +- **Misleading warnings at the end:** S1Tiling may report "download failures" + for redundant S1 products that weren't strictly needed. Check if the actual + output files were produced. + +--- + +## 7 – Expected output files + +After a successful run, you should find: + +### Final products: `$S1T_WORKDIR/data_out/31TCH/` + +| Pattern | Description | +|---|---| +| `s1a_31TCH_vv_DES_{orbit}_{date}_GammaNaughtRTC.tif` | γ0T VV backscatter | +| `s1a_31TCH_vh_DES_{orbit}_{date}_GammaNaughtRTC.tif` | γ0T VH backscatter | +| `s1a_31TCH_vv_DES_{orbit}_{date}_GammaNaughtRTC_BorderMask.tif` | VV border mask | +| `s1a_31TCH_vh_DES_{orbit}_{date}_GammaNaughtRTC_BorderMask.tif` | VH border mask | + +### Gamma Area maps: `$S1T_WORKDIR/data_gamma_area/` + +| Pattern | Description | +|---|---| +| `GAMMA_AREA_s1a_31TCH_DES_{orbit}.tif` | γ area map for this S2 tile + orbit | + +### GeoTIFF metadata (from S1Tiling docs) + +Each final product should contain these GeoTIFF tags: + +| Tag | Example value | +|---|---| +| `CALIBRATION` | `GammaNaughtRTC` | +| `IMAGE_TYPE` | `BACKSCATTERING` | +| `FLYING_UNIT_CODE` | `s1a` | +| `POLARIZATION` | `vv` or `vh` | +| `ORBIT_DIRECTION` | `DES` | +| `RELATIVE_ORBIT_NUMBER` | e.g. `110` | +| `S2_TILE_CORRESPONDING_CODE` | `31TCH` | +| `SPATIAL_RESOLUTION` | `10.0` | +| `ORTHORECTIFIED` | `true` | +| `GAMMA_AREA_FILE` | name of the gamma area map used | +| `TIFFTAG_SOFTWARE` | `S1 Tiling v1.4.0` | +| `ACQUISITION_DATETIME` | UTC datetime | + +Product encoding: **Float32 GeoTIFF, deflate compressed**, CRS matching the +MGRS tile UTM zone (e.g., EPSG:32631 for 31TCH). + +--- + +## 8 – Inspect the outputs + +Run the inspector script on the results: + +```bash +cd /home/emathot/Workspace/eopf-explorer/data-model + +# Inspect all final products with validation +python analysis/inspect_s1tiling_geotiff.py --validate "$S1T_WORKDIR/data_out/" + +# Inspect gamma area maps +python analysis/inspect_s1tiling_geotiff.py --validate "$S1T_WORKDIR/data_gamma_area/" + +# JSON output for programmatic analysis +python analysis/inspect_s1tiling_geotiff.py --json "$S1T_WORKDIR/data_out/" > analysis/s1tiling_output_metadata.json +``` + +### What to report back + +Please share the following after the run: + +1. **Console output** of the docker run (especially the final execution report) +2. **File listing:** + ```bash + find "$S1T_WORKDIR"/data_out -name "*.tif" -ls + find "$S1T_WORKDIR"/data_gamma_area -name "*.tif" -ls + ``` +3. **Inspector output:** + ```bash + python analysis/inspect_s1tiling_geotiff.py --validate "$S1T_WORKDIR/data_out/" + python analysis/inspect_s1tiling_geotiff.py --validate "$S1T_WORKDIR/data_gamma_area/" + ``` +4. **JSON dump** (for detailed programmatic review): + ```bash + python analysis/inspect_s1tiling_geotiff.py --json "$S1T_WORKDIR/data_out/" \ + "$S1T_WORKDIR/data_gamma_area/" > analysis/s1tiling_output_metadata.json + ``` + +--- + +## 9 – Alternative: different MGRS tile / shorter run + +If 31TCH doesn't work (DEM unavailable, no S1 coverage in the window, etc.), +pick another tile. Good candidates with frequent S1 coverage: + +| MGRS Tile | Location | UTM Zone | +|---|---|---| +| `31TCH` | South France (Toulouse area) | 31N | +| `32TQM` | North Italy | 32N | +| `33UUP` | Germany | 33N | +| `10SEG` | California coast | 10N | + +Adjust `tiles`, `roi_by_tiles`, and DEM tiles accordingly. + +To further reduce processing time, try `calibration: gamma` (simple γ0 without +RTC). This skips the Gamma Area map computation but won't produce the RTC +product we actually need. + +--- + +## 10 – Filename pattern for the GAMMA_AREA file + +Note from the docs: the Gamma Area filename uses the format +`GAMMA_AREA_s1a_{tile}_{direction}_{orbit}.tif` (example: +`GAMMA_AREA_s1a_31TCH_DES_110.tif`). This is different from the documented +template `GAMMA_AREA_{tile}_{orbit}.tif` — the actual filename includes the +flying unit code and orbit direction. The inspector script handles both +patterns. diff --git a/analysis/s1tiling_eodag4_patch.py b/analysis/s1tiling_eodag4_patch.py new file mode 100644 index 00000000..e9a3e4ca --- /dev/null +++ b/analysis/s1tiling_eodag4_patch.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python3 +""" +Patch S1Tiling 1.4.0 for EODAG 4.0.0 compatibility. + +EODAG 4.0.0 introduced several breaking changes for S1Tiling: +1. `productType` kwarg to dag.search() was renamed to `collection`. + Having both causes cop_dataspace to fail silently → falls back to peps. +2. Product properties use STAC names (sat:orbit_state, platform, etc.) + instead of legacy EODAG names (orbitDirection, platformSerialIdentifier, etc.) +3. cop_dataspace OData v4 rejects `polarizationChannels` and `sensorMode`. +4. cop_dataspace requires UPPERCASE orbit direction ("DESCENDING" not "descending"). +5. `relativeOrbitNumber` search param silently returns 0 results on cop_dataspace. + +This script patches: +- S1FileManager.py: fixes the search() call (issues 1, 3, 4, 5) +- s1/product.py: adds legacy→STAC property name fallback (issue 2) + +Usage (inside the Docker container): + python3 /patch/s1tiling_eodag4_patch.py +""" + +import pathlib +import re + +S1T_PKG = pathlib.Path("/opt/S1TilingEnv/lib/python3.10/site-packages/s1tiling/libs") + + +def patch_s1filemanager(): + """Fix search() call: add collection param, remove unsupported kwargs.""" + fpath = S1T_PKG / "S1FileManager.py" + src = fpath.read_text() + + # Replace productType with collection (EODAG 4.0.0 rename). + # Keeping productType alongside collection causes cop_dataspace to fail + # silently, making EODAG fall back to peps. + src = src.replace( + "productType=product_type,", + "collection=product_type,", + 1, # only first occurrence + ) + + # Remove polarizationChannels (unsupported by cop_dataspace OData) + src = re.sub(r"\n\s*# If we have eodag.*\n", "\n", src) + src = re.sub(r"\n\s*polarizationChannels=dag_polarization_param,", "", src) + + # Remove sensorMode="IW" (unsupported by cop_dataspace OData) + src = re.sub(r'\n\s*sensorMode="IW",', "", src) + + # Remove relativeOrbitNumber (not supported by cop_dataspace OData — + # returns 0 results silently). S1Tiling has post-search filtering for + # relative orbits when len(relative_orbit_list) > 1, and when list + # has exactly 1 element, orbitNumber=None was passed anyway for most configs. + src = re.sub( + r"\n\s*relativeOrbitNumber=dag_orbit_list_param,.*", + "", + src, + ) + + # cop_dataspace OData requires UPPERCASE orbit direction values + # ("DESCENDING" not "descending"). S1Tiling's k_dir_assoc produces lowercase. + src = src.replace( + "{ 'ASC': 'ascending', 'DES': 'descending' }", + "{ 'ASC': 'ASCENDING', 'DES': 'DESCENDING' }", + ) + + fpath.write_text(src) + print(f" Patched {fpath.name}") + + +def patch_product_property(): + """Add STAC property name fallback to product_property().""" + fpath = S1T_PKG / "s1" / "product.py" + src = fpath.read_text() + + old = ( + "def product_property(prod: EOProduct, key: str, default=None):\n" + ' """\n' + " Returns the required (EODAG) product property, " + "or default in the property isn't found.\n" + ' """\n' + " res = prod.properties.get(key, default)\n" + " return res" + ) + + new = ( + "def product_property(prod: EOProduct, key: str, default=None):\n" + ' """\n' + " Returns the required (EODAG) product property, " + "or default in the property isn't found.\n" + " EODAG 4.0.0 uses STAC property names; " + "fall back to them for legacy keys.\n" + ' """\n' + " _FALLBACK = {\n" + ' "orbitDirection": "sat:orbit_state",\n' + ' "platformSerialIdentifier": "platform",\n' + ' "relativeOrbitNumber": "sat:relative_orbit",\n' + ' "orbitNumber": "sat:absolute_orbit",\n' + ' "polarizationChannels": "sar:polarizations",\n' + ' "startTimeFromAscendingNode": "start_datetime",\n' + ' "completionTimeFromAscendingNode": "end_datetime",\n' + " }\n" + " res = prod.properties.get(key, None)\n" + " if res is None and key in _FALLBACK:\n" + " res = prod.properties.get(_FALLBACK[key], default)\n" + ' if key == "polarizationChannels" and isinstance(res, list):\n' + ' res = "+".join(res)\n' + " return res if res is not None else default" + ) + + if old not in src: + print(f" WARNING: product_property() not found in {fpath.name}, skipping") + return + + src = src.replace(old, new) + fpath.write_text(src) + print(f" Patched {fpath.name}") + + +if __name__ == "__main__": + print("Applying S1Tiling EODAG 4.0.0 compatibility patches...") + patch_s1filemanager() + patch_product_property() + print("Done.") diff --git a/analysis/s1tiling_output_metadata.json b/analysis/s1tiling_output_metadata.json new file mode 100644 index 00000000..6ca5266b --- /dev/null +++ b/analysis/s1tiling_output_metadata.json @@ -0,0 +1,998 @@ +[ + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vh_DES_008_20250210t060920_GammaNaughtRTC.tif", + "filename": "s1a_31TCH_vh_DES_008_20250210t060920_GammaNaughtRTC.tif", + "file_type": "gamma_naught_rtc", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vh", + "orbit_direction": "DES", + "orbit": "008", + "acquisition_stamp": "20250210t060920" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "float32", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Gamma0 RTC Calibrated Sentinel-1A IW GRD", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:31:27", + "ACQUISITION_DATETIME": "2025:02:10T06:09:20Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_008.tif", + "IMAGE_TYPE": "GRD", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250210T060920_20250210T060945_057830_0721D8_1F28", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057830", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vh", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-10T06:50:29.028854Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "008", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26497", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vh_DES_008_20250210t060920_GammaNaughtRTC_BorderMask.tif", + "filename": "s1a_31TCH_vh_DES_008_20250210t060920_GammaNaughtRTC_BorderMask.tif", + "file_type": "border_mask", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vh", + "orbit_direction": "DES", + "orbit": "008", + "acquisition_stamp": "20250210t060920" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "uint8", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Orthorectified Sentinel-1A IW GRD smoothed border mask S2 tile", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:31:37", + "ACQUISITION_DATETIME": "2025:02:10T06:09:20Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_008.tif", + "IMAGE_TYPE": "MASK", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250210T060920_20250210T060945_057830_0721D8_1F28", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057830", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vh", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-10T06:50:29.028854Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "008", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26497", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vh_DES_037_20250212t055301_GammaNaughtRTC.tif", + "filename": "s1a_31TCH_vh_DES_037_20250212t055301_GammaNaughtRTC.tif", + "file_type": "gamma_naught_rtc", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vh", + "orbit_direction": "DES", + "orbit": "037", + "acquisition_stamp": "20250212t055301" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "float32", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Gamma0 RTC Calibrated Sentinel-1A IW GRD", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:34:20", + "ACQUISITION_DATETIME": "2025:02:12T05:53:01Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_037.tif", + "IMAGE_TYPE": "GRD", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250212T055301_20250212T055326_057859_0722FE_13B0", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057859", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vh", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-12T06:24:23.751945Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "037", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26399", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vh_DES_037_20250212t055301_GammaNaughtRTC_BorderMask.tif", + "filename": "s1a_31TCH_vh_DES_037_20250212t055301_GammaNaughtRTC_BorderMask.tif", + "file_type": "border_mask", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vh", + "orbit_direction": "DES", + "orbit": "037", + "acquisition_stamp": "20250212t055301" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "uint8", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Orthorectified Sentinel-1A IW GRD smoothed border mask S2 tile", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:34:30", + "ACQUISITION_DATETIME": "2025:02:12T05:53:01Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_037.tif", + "IMAGE_TYPE": "MASK", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250212T055301_20250212T055326_057859_0722FE_13B0", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057859", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vh", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-12T06:24:23.751945Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "037", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26399", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vh_DES_110_20250205t060110_GammaNaughtRTC.tif", + "filename": "s1a_31TCH_vh_DES_110_20250205t060110_GammaNaughtRTC.tif", + "file_type": "gamma_naught_rtc", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vh", + "orbit_direction": "DES", + "orbit": "110", + "acquisition_stamp": "20250205t060110" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "float32", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Gamma0 RTC Calibrated Sentinel-1A IW GRD", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:38:03", + "ACQUISITION_DATETIME": "2025:02:05T06:01:10Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_110.tif", + "IMAGE_TYPE": "GRD", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250205T060110_20250205T060135_057757_071EE7_D5AA", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057757", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vh", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-05T06:39:31.395792Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "110", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26416", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vh_DES_110_20250205t060110_GammaNaughtRTC_BorderMask.tif", + "filename": "s1a_31TCH_vh_DES_110_20250205t060110_GammaNaughtRTC_BorderMask.tif", + "file_type": "border_mask", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vh", + "orbit_direction": "DES", + "orbit": "110", + "acquisition_stamp": "20250205t060110" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "uint8", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Orthorectified Sentinel-1A IW GRD smoothed border mask S2 tile", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:38:14", + "ACQUISITION_DATETIME": "2025:02:05T06:01:10Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_110.tif", + "IMAGE_TYPE": "MASK", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250205T060110_20250205T060135_057757_071EE7_D5AA", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057757", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vh", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-05T06:39:31.395792Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "110", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26416", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vv_DES_008_20250210t060920_GammaNaughtRTC.tif", + "filename": "s1a_31TCH_vv_DES_008_20250210t060920_GammaNaughtRTC.tif", + "file_type": "gamma_naught_rtc", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vv", + "orbit_direction": "DES", + "orbit": "008", + "acquisition_stamp": "20250210t060920" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "float32", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Gamma0 RTC Calibrated Sentinel-1A IW GRD", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:31:25", + "ACQUISITION_DATETIME": "2025:02:10T06:09:20Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_008.tif", + "IMAGE_TYPE": "GRD", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250210T060920_20250210T060945_057830_0721D8_1F28", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057830", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vv", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-10T06:50:29.028854Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "008", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26497", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vv_DES_008_20250210t060920_GammaNaughtRTC_BorderMask.tif", + "filename": "s1a_31TCH_vv_DES_008_20250210t060920_GammaNaughtRTC_BorderMask.tif", + "file_type": "border_mask", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vv", + "orbit_direction": "DES", + "orbit": "008", + "acquisition_stamp": "20250210t060920" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "uint8", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Orthorectified Sentinel-1A IW GRD smoothed border mask S2 tile", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:31:29", + "ACQUISITION_DATETIME": "2025:02:10T06:09:20Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_008.tif", + "IMAGE_TYPE": "MASK", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250210T060920_20250210T060945_057830_0721D8_1F28", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057830", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vv", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-10T06:50:29.028854Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "008", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26497", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vv_DES_037_20250212t055301_GammaNaughtRTC.tif", + "filename": "s1a_31TCH_vv_DES_037_20250212t055301_GammaNaughtRTC.tif", + "file_type": "gamma_naught_rtc", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vv", + "orbit_direction": "DES", + "orbit": "037", + "acquisition_stamp": "20250212t055301" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "float32", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Gamma0 RTC Calibrated Sentinel-1A IW GRD", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:34:19", + "ACQUISITION_DATETIME": "2025:02:12T05:53:01Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_037.tif", + "IMAGE_TYPE": "GRD", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250212T055301_20250212T055326_057859_0722FE_13B0", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057859", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vv", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-12T06:24:23.751945Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "037", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26399", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vv_DES_037_20250212t055301_GammaNaughtRTC_BorderMask.tif", + "filename": "s1a_31TCH_vv_DES_037_20250212t055301_GammaNaughtRTC_BorderMask.tif", + "file_type": "border_mask", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vv", + "orbit_direction": "DES", + "orbit": "037", + "acquisition_stamp": "20250212t055301" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "uint8", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Orthorectified Sentinel-1A IW GRD smoothed border mask S2 tile", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:34:22", + "ACQUISITION_DATETIME": "2025:02:12T05:53:01Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_037.tif", + "IMAGE_TYPE": "MASK", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250212T055301_20250212T055326_057859_0722FE_13B0", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057859", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vv", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-12T06:24:23.751945Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "037", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26399", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vv_DES_110_20250205t060110_GammaNaughtRTC.tif", + "filename": "s1a_31TCH_vv_DES_110_20250205t060110_GammaNaughtRTC.tif", + "file_type": "gamma_naught_rtc", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vv", + "orbit_direction": "DES", + "orbit": "110", + "acquisition_stamp": "20250205t060110" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "float32", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Gamma0 RTC Calibrated Sentinel-1A IW GRD", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:38:01", + "ACQUISITION_DATETIME": "2025:02:05T06:01:10Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_110.tif", + "IMAGE_TYPE": "GRD", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250205T060110_20250205T060135_057757_071EE7_D5AA", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057757", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vv", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-05T06:39:31.395792Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "110", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26416", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + }, + { + "file": "/home/emathot/Downloads/s1tiling_test/data_out/31TCH/s1a_31TCH_vv_DES_110_20250205t060110_GammaNaughtRTC_BorderMask.tif", + "filename": "s1a_31TCH_vv_DES_110_20250205t060110_GammaNaughtRTC_BorderMask.tif", + "file_type": "border_mask", + "filename_metadata": { + "flying_unit_code": "s1a", + "tile_name": "31TCH", + "polarisation": "vv", + "orbit_direction": "DES", + "orbit": "110", + "acquisition_stamp": "20250205t060110" + }, + "rasterio_profile": { + "driver": "GTiff", + "dtype": "uint8", + "width": 10980, + "height": 10980, + "count": 1, + "crs": "EPSG:32631", + "nodata": 0.0 + }, + "transform": { + "affine": [ + 10.0, + 0.0, + 299999.9999974121, + 0.0, + -10.0, + 4799999.99999915 + ], + "pixel_size_x": 10.0, + "pixel_size_y": 10.0 + }, + "bounds": { + "left": 299999.9999974121, + "bottom": 4690199.99999915, + "right": 409799.9999974121, + "top": 4799999.99999915 + }, + "tags": { + "TIFFTAG_IMAGEDESCRIPTION": "Orthorectified Sentinel-1A IW GRD smoothed border mask S2 tile", + "TIFFTAG_SOFTWARE": "S1 Tiling v1.4.0", + "TIFFTAG_DATETIME": "2026:03:23 06:38:06", + "ACQUISITION_DATETIME": "2025:02:05T06:01:10Z", + "CALIBRATION": "GammaNaughtRTC", + "DataType": "3", + "DEM_INFO": "SRTM 30m", + "FACILITY_IDENTIFIER": "ESA Sentinel-1 IPF 003.90", + "FLYING_UNIT_CODE": "s1a", + "GAMMA_AREA_FILE": "GAMMA_AREA_31TCH_110.tif", + "IMAGE_TYPE": "MASK", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20250205T060110_20250205T060135_057757_071EE7_D5AA", + "LOWER_SIGNAL_VALUE": "1e-07", + "METADATATYPE": "OTB", + "NoData": "0", + "NOISE_REMOVED": "True", + "ORBIT_DIRECTION": "DES", + "ORBIT_NUMBER": "057757", + "ORTHORECTIFICATION_INTERPOLATOR": "bco", + "ORTHORECTIFIED": "true", + "OTB_VERSION": "9.1.1", + "phase0_test": "s1-grd-rtc-validation", + "POLARIZATION": "vv", + "producer": "eopf-geozarr-phase0", + "ProductionDate": "2025-02-05T06:39:31.395792Z", + "ProductType": "GRD", + "RELATIVE_ORBIT_NUMBER": "110", + "S2_TILE_CORRESPONDING_CODE": "31TCH", + "SPATIAL_RESOLUTION": "10.0", + "TileHintX": "26416", + "TileHintY": "1", + "AREA_OR_POINT": "Area" + }, + "band_descriptions": [ + null + ], + "band_tags": { + "band_1": { + "NoData": "0", + "BandName": "\u03b3\u00b0RTC" + } + } + } +] diff --git a/issues/s1tiling-eodag-collection-bug.md b/issues/s1tiling-eodag-collection-bug.md new file mode 100644 index 00000000..5ba45ea9 --- /dev/null +++ b/issues/s1tiling-eodag-collection-bug.md @@ -0,0 +1,142 @@ +# S1Tiling EODAG 4.0 compatibility: `productType` deprecated, `collection` now required + +## Title + +S1 product search fails with EODAG ≥ 4.0: `productType` is deprecated, `collection` parameter required + +## Description + +Starting with EODAG 4.0.0, the `productType` keyword argument to `dag.search()` was replaced by `collection`. The old `productType` parameter is no longer recognized and `collection` is now mandatory. This affects S1Tiling 1.4.0 (Docker image `registry.orfeo-toolbox.org/s1-tiling/s1tiling:1.4.0-ubuntu-otb9.1.1`) which ships EODAG 4.0.0, and the bug persists on the current `develop` branch (version 2.0.0). + +In `s1tiling/libs/s1/file_manager.py`, the `as_eodag_parameters()` method (line 195) builds a search dict with the key `"productType"`. This dict is unpacked via `**search_args` in the `dag.search()` call (line 443). EODAG 4.0.0 raises `ValidationError("Field required: collection")` because the deprecated `productType` kwarg is not recognized and the required `collection` parameter is missing. + +> **Note:** The orbit provider (`s1tiling/libs/orbit/_providers.py`, line 199–202) already uses `collection="SENTINEL-1"` alongside `productType="S1_AUX_POEORB"` in `search_all()`, so the orbit file search is already partially adapted. Only the main S1 product search is broken. + +### Verified fix + +Tested inside the 1.4.0 Docker container: +```python +from eodag import EODataAccessGateway +dag = EODataAccessGateway() + +# FAILS — productType no longer recognized: +dag.search(productType="S1_SAR_GRD", limit=1, box=(0, 42, 2, 44)) +# → ValidationError: Field required: collection + +# WORKS — collection is the new parameter name: +dag.search(collection="S1_SAR_GRD", limit=1, box=(0, 42, 2, 44)) +# → 1 result from cop_dataspace +``` + +## Steps to Reproduce + +1. Pull the 1.4.0 Docker image: + ```bash + docker pull registry.orfeo-toolbox.org/s1-tiling/s1tiling:1.4.0-ubuntu-otb9.1.1 + ``` + +2. Configure `eodag.yml` with valid `cop_dataspace` credentials. + +3. Run `S1Processor` with `download: True` and any tile configuration: + ``` + WARNING - Cannot download S1 images associated to 31TCH: + Cannot request products for tile 31TCH on data provider: Field required: collection + ``` + +## Root Cause + +In `s1tiling/libs/s1/file_manager.py`, the `SearchCriteria.as_eodag_parameters()` method returns: + +```python +# s1tiling/libs/s1/file_manager.py, line 193–205 +res = { + "productType" : self.__product_type, # ← deprecated in EODAG 4.0 + "start" : self.__first_date, + "end" : self.__last_date, + "sensorMode" : self.__sensor_mode, + "polarizationChannels" : dag_polarization_param, + "orbitDirection" : dag_orbit_dir_param, + "relativeOrbitNumber" : dag_orbit_list_param, + "platformSerialIdentifier": dag_platform_list_param, +} +``` + +This dict is then unpacked into the search call: + +```python +# s1tiling/libs/s1/file_manager.py, line 443–446 +page_products = dag.search( + page=page, items_per_page=self.__searched_items_per_page, + raise_errors=True, + geom=footprint, + **search_args, # ← includes productType, sensorMode, polarizationChannels +) +``` + +EODAG 4.0 renamed the `productType` parameter to `collection` (the value `"S1_SAR_GRD"` stays the same). + +## Suggested Fix + +Two issues need fixing in `as_eodag_parameters()`: + +### 1. Replace `"productType"` key with `"collection"` + +In `s1tiling/libs/s1/file_manager.py`, line 195, change: +```python +"productType" : self.__product_type, +``` +to: +```python +"collection" : self.__product_type, +``` + +### 2. Remove `"polarizationChannels"` and `"sensorMode"` entries + +When using `cop_dataspace`, the OData v4 API rejects `polarizationChannels` and `sensorMode` as invalid fields (HTTP 400: `"Invalid field: polarizationChannels"`). These parameters are not mapped in EODAG's `cop_dataspace` provider config and are passed through raw to the OData endpoint, which rejects them. + +`orbitDirection`, `relativeOrbitNumber`, and `platformSerialIdentifier` are fine — they have proper mappings in EODAG's provider config. + +The `filter_eodag_search_results()` method already handles post-search filtering for polarization and platform, but `sensorMode` is currently only enforced at the search level. It should be moved to post-search filtering as well (or guarded conditionally per provider). + +Remove both entries from the returned dict: +```python +res = { + "collection" : self.__product_type, + "start" : self.__first_date, + "end" : self.__last_date, + # "sensorMode" : self.__sensor_mode, # removed: rejected by cop_dataspace OData + # "polarizationChannels" : dag_polarization_param, # removed: rejected by cop_dataspace OData + "orbitDirection" : dag_orbit_dir_param, + "relativeOrbitNumber" : dag_orbit_list_param, + "platformSerialIdentifier": dag_platform_list_param, +} +``` + +> **Note:** Removing `sensorMode` from the search parameters means the S1 query will no longer filter by acquisition mode server-side. Since S1Tiling only targets IW-mode GRD products and `productType="S1_SAR_GRD"` already constrains results to GRD, this is unlikely to cause issues in practice, but a post-search filter on `sensorMode` should be added to `filter_eodag_search_results()` for correctness. + +## Workaround + +Patch the file inside the Docker container before running: + +```bash +docker exec -it python3 -c " +import pathlib +p = pathlib.Path('/opt/S1TilingEnv/lib/python3.10/site-packages/s1tiling/libs/s1/file_manager.py') +s = p.read_text() +s = s.replace('\"productType\"', '\"collection\"', 1) +s = s.replace('\"sensorMode\" : self.__sensor_mode,\n', '') +s = s.replace('\"polarizationChannels\" : dag_polarization_param,\n', '') +p.write_text(s) +" +``` + +## Environment + +- **S1Tiling:** 1.4.0 (Docker) / 2.0.0 (`develop` — bug persists) +- **EODAG:** ≥ 4.0.0 +- **Docker image:** `registry.orfeo-toolbox.org/s1-tiling/s1tiling:1.4.0-ubuntu-otb9.1.1` +- **Provider:** `cop_dataspace` (Copernicus Data Space Ecosystem) + +## Labels + +`bug`, `eodag` diff --git a/pyproject.toml b/pyproject.toml index 2e2ce41d..47a8017e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -149,6 +149,9 @@ ignore = [ "TRY003", # Long exception messages outside class - common pattern ] +[tool.ruff.lint.flake8-type-checking] +runtime-evaluated-base-classes = ["pydantic.BaseModel"] + [tool.pyright] include = ["src", "tests"] pythonVersion = "3.12" diff --git a/src/eopf_geozarr/cli.py b/src/eopf_geozarr/cli.py index 5e4d95da..ea4bc867 100755 --- a/src/eopf_geozarr/cli.py +++ b/src/eopf_geozarr/cli.py @@ -1054,6 +1054,60 @@ def validate_command(args: argparse.Namespace) -> None: sys.exit(1) +# ============================================================================= +# S1 Ingestion Commands +# ============================================================================= + + +def ingest_s1_command(args: argparse.Namespace) -> None: + """Ingest a single S1Tiling acquisition into a GeoZarr V3 store.""" + from .conversion.s1_ingest import ingest_s1tiling_acquisition + + try: + idx = ingest_s1tiling_acquisition( + vv_path=args.vv, + vh_path=args.vh, + border_mask_path=args.mask, + store_path=args.store, + orbit_direction=args.orbit_dir, + ) + log.info("✅ Acquisition ingested", time_index=idx, store=args.store) + except Exception as e: + log.error("❌ Error ingesting acquisition", error=str(e)) + sys.exit(1) + + +def ingest_s1_conditions_command(args: argparse.Namespace) -> None: + """Ingest S1Tiling condition arrays into a GeoZarr V3 store.""" + from .conversion.s1_ingest import ingest_s1tiling_conditions + + try: + ingest_s1tiling_conditions( + store_path=args.store, + orbit_direction=args.orbit_dir, + relative_orbit=args.relative_orbit, + gamma_area_path=getattr(args, "gamma_area", None), + lia_path=getattr(args, "lia", None), + incidence_angle_path=getattr(args, "incidence_angle", None), + ) + log.info("✅ Conditions ingested", store=args.store, orbit_dir=args.orbit_dir) + except Exception as e: + log.error("❌ Error ingesting conditions", error=str(e)) + sys.exit(1) + + +def consolidate_s1_command(args: argparse.Namespace) -> None: + """Consolidate metadata for an S1 GeoZarr store.""" + from .conversion.s1_ingest import consolidate_s1_store + + try: + consolidate_s1_store(args.store, args.orbit_dir) + log.info("✅ Metadata consolidated", store=args.store) + except Exception as e: + log.error("❌ Error consolidating metadata", error=str(e)) + sys.exit(1) + + def create_parser() -> argparse.ArgumentParser: """ Create the argument parser for the CLI. @@ -1182,9 +1236,74 @@ def create_parser() -> argparse.ArgumentParser: # Add S2 optimization commands add_s2_optimization_commands(subparsers) + # Add S1 ingestion commands + add_s1_ingestion_commands(subparsers) + return parser +def add_s1_ingestion_commands(subparsers: argparse._SubParsersAction) -> None: + """Add S1 GRD RTC ingestion commands to CLI parser.""" + + # ingest-s1: single acquisition + s1_parser = subparsers.add_parser( + "ingest-s1", help="Ingest a single S1Tiling acquisition into a GeoZarr V3 store" + ) + s1_parser.add_argument("--vv", type=str, required=True, help="Path to VV GeoTIFF") + s1_parser.add_argument("--vh", type=str, required=True, help="Path to VH GeoTIFF") + s1_parser.add_argument("--mask", type=str, required=True, help="Path to border mask GeoTIFF") + s1_parser.add_argument("--store", type=str, required=True, help="Path to output Zarr V3 store") + s1_parser.add_argument( + "--orbit-dir", + type=str, + required=True, + choices=["ascending", "descending"], + help="Orbit direction", + ) + s1_parser.set_defaults(func=ingest_s1_command) + + # ingest-s1-conditions: condition arrays + cond_parser = subparsers.add_parser( + "ingest-s1-conditions", + help="Ingest S1Tiling condition arrays (gamma_area, LIA) into a GeoZarr V3 store", + ) + cond_parser.add_argument( + "--store", type=str, required=True, help="Path to existing Zarr V3 store" + ) + cond_parser.add_argument( + "--orbit-dir", + type=str, + required=True, + choices=["ascending", "descending"], + help="Orbit direction", + ) + cond_parser.add_argument( + "--relative-orbit", type=int, required=True, help="Relative orbit number (e.g. 37)" + ) + cond_parser.add_argument( + "--gamma-area", type=str, default=None, help="Path to gamma area GeoTIFF" + ) + cond_parser.add_argument("--lia", type=str, default=None, help="Path to LIA GeoTIFF") + cond_parser.add_argument( + "--incidence-angle", type=str, default=None, help="Path to incidence angle GeoTIFF" + ) + cond_parser.set_defaults(func=ingest_s1_conditions_command) + + # consolidate-s1: metadata consolidation + cons_parser = subparsers.add_parser( + "consolidate-s1", help="Consolidate metadata for an S1 GeoZarr V3 store" + ) + cons_parser.add_argument("--store", type=str, required=True, help="Path to Zarr V3 store") + cons_parser.add_argument( + "--orbit-dir", + type=str, + required=True, + choices=["ascending", "descending"], + help="Orbit direction", + ) + cons_parser.set_defaults(func=consolidate_s1_command) + + def add_s2_optimization_commands(subparsers: argparse._SubParsersAction) -> None: """Add S2 optimization commands to CLI parser.""" diff --git a/src/eopf_geozarr/conversion/__init__.py b/src/eopf_geozarr/conversion/__init__.py index 1a0e448a..b58b2c54 100644 --- a/src/eopf_geozarr/conversion/__init__.py +++ b/src/eopf_geozarr/conversion/__init__.py @@ -17,6 +17,14 @@ iterative_copy, setup_datatree_metadata_geozarr_spec_compliant, ) +from .s1_ingest import ( + consolidate_s1_store, + discover_s1tiling_acquisitions, + discover_s1tiling_conditions, + extract_geotiff_metadata, + ingest_s1tiling_acquisition, + ingest_s1tiling_conditions, +) from .utils import ( calculate_aligned_chunk_size, downsample_2d_array, @@ -29,9 +37,15 @@ "calculate_aligned_chunk_size", "calculate_overview_levels", "consolidate_metadata", + "consolidate_s1_store", "create_geozarr_dataset", + "discover_s1tiling_acquisitions", + "discover_s1tiling_conditions", "downsample_2d_array", + "extract_geotiff_metadata", "get_s3_credentials_info", + "ingest_s1tiling_acquisition", + "ingest_s1tiling_conditions", "is_grid_mapping_variable", "is_s3_path", "iterative_copy", diff --git a/src/eopf_geozarr/conversion/s1_ingest.py b/src/eopf_geozarr/conversion/s1_ingest.py new file mode 100644 index 00000000..8109e8a7 --- /dev/null +++ b/src/eopf_geozarr/conversion/s1_ingest.py @@ -0,0 +1,930 @@ +"""S1 GRD RTC GeoTIFF → GeoZarr V3 ingestion pipeline. + +Converts S1Tiling γ0T RTC GeoTIFF outputs into a sharded Zarr V3 store +with multiscale overviews, spatial coordinate arrays, and full GeoZarr +convention metadata. + +Public API: + - extract_geotiff_metadata(path) -> S1TilingMetadata + - ingest_s1tiling_acquisition(vv_path, vh_path, border_mask_path, store_path, orbit_direction) -> int + - ingest_s1tiling_conditions(store_path, orbit_direction, relative_orbit, ...) -> None + - consolidate_s1_store(store_path, orbit_direction) -> None + - discover_s1tiling_acquisitions(input_dir) -> list[dict] + - discover_s1tiling_conditions(input_dir) -> list[dict] +""" + +from __future__ import annotations + +import re +from dataclasses import dataclass +from math import ceil +from pathlib import Path + +import numpy as np +import rasterio +import structlog +import zarr +import zarr.codecs +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 + +log = structlog.get_logger() + +# ============================================================================= +# Constants +# ============================================================================= + +MULTISCALES_UUID = multiscales_cm.UUID +GEO_PROJ_UUID = geo_proj.UUID +SPATIAL_UUID = spatial_cm.UUID + +ZARR_CONVENTIONS = [multiscales_cm.CMO, geo_proj.CMO, spatial_cm.CMO] + +# Overview chain: (level_name, parent_name, downsample_factor) +OVERVIEW_CHAIN = [ + ("r10m", None, 1), + ("r20m", "r10m", 2), + ("r60m", "r20m", 3), + ("r120m", "r60m", 2), + ("r360m", "r120m", 3), + ("r720m", "r360m", 2), +] + +# S1Tiling filename pattern +# e.g. s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC.tif +# Multi-frame products mask the time as 'txxxxxx' (no single shared acquisition time); those are +# accepted here and the real stamp is resolved from the GeoTIFF ACQUISITION_DATETIME tag (#183). +S1TILING_FILENAME_PATTERN = re.compile( + r"(?Ps1[abc])_" + r"(?P[0-9]{2}[A-Z]{3})_" + r"(?Pvv|vh)_" + r"(?PASC|DES)_" + r"(?P\d{3})_" + r"(?P\d{8}t(?:\d{6}|x{6}))_" + r"(?PGammaNaughtRTC)" + r"(?P_BorderMask)?\.tif$" +) + +# S1Tiling conditions filename patterns +# e.g. GAMMA_AREA_31TCH_008.tif or GAMMA_AREA_s1a_31TCH_ASC_008.tif +S1TILING_GAMMA_AREA_PATTERN = re.compile( + r"^GAMMA_AREA_(?:s1[abc]_)?(?P[A-Z0-9]+)_(?:(?:ASC|DES)_)?(?P\d{3})\.tif$", + re.IGNORECASE, +) +# e.g. sin_LIA_31TCH_008.tif or LIA_31TCH_008.tif +S1TILING_LIA_PATTERN = re.compile( + r"^(?Psin_LIA|LIA)_(?P[A-Z0-9]+)_(?P\d{3})\.tif$", + re.IGNORECASE, +) + + +# ============================================================================= +# Data Transfer Object +# ============================================================================= + + +@dataclass(frozen=True) +class S1TilingMetadata: + """Metadata extracted from an S1Tiling GeoTIFF.""" + + crs: str + spatial_transform: list[float] + shape: list[int] + bounds: list[float] + datetime: str + absolute_orbit: int + relative_orbit: int + platform: str + calibration: str + input_s1_images: str + + +# ============================================================================= +# Metadata Extraction +# ============================================================================= + + +def _normalise_s1tiling_datetime(dt_str: str) -> str: + """Normalise S1Tiling datetime format to ISO 8601. + + Input: "2025:02:10T06:09:20Z" (S1Tiling uses colons in date part) + Output: "2025-02-10T06:09:20" + """ + dt_normalised = dt_str.replace("Z", "") + parts = dt_normalised.split("T") + if len(parts) == 2: + date_part = parts[0].replace(":", "-") + dt_normalised = f"{date_part}T{parts[1]}" + return dt_normalised + + +def extract_geotiff_metadata(path: str | Path) -> S1TilingMetadata: + """Extract CRS, transform, bounds, and custom tags from an S1Tiling GeoTIFF. + + Raises + ------ + ValueError + If critical tags (ACQUISITION_DATETIME, ORBIT_NUMBER, + RELATIVE_ORBIT_NUMBER, FLYING_UNIT_CODE) are missing. + """ + with 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] + + # Validate critical tags + required_tags = [ + "ACQUISITION_DATETIME", + "ORBIT_NUMBER", + "RELATIVE_ORBIT_NUMBER", + "FLYING_UNIT_CODE", + ] + missing = [tag for tag in required_tags if tag not in tags] + if missing: + raise ValueError(f"GeoTIFF {path} missing required tags: {missing}") + + dt_raw = tags["ACQUISITION_DATETIME"] + dt_normalised = _normalise_s1tiling_datetime(dt_raw) + + metadata = S1TilingMetadata( + crs=str(src.crs), + spatial_transform=spatial_transform, + shape=[src.height, src.width], + bounds=[src.bounds.left, src.bounds.bottom, src.bounds.right, src.bounds.top], + datetime=dt_normalised, + absolute_orbit=int(tags["ORBIT_NUMBER"]), + relative_orbit=int(tags["RELATIVE_ORBIT_NUMBER"]), + platform=tags["FLYING_UNIT_CODE"], + calibration=tags.get("CALIBRATION", ""), + input_s1_images=tags.get("INPUT_S1_IMAGES", ""), + ) + + log.info( + "Extracted GeoTIFF metadata", + path=str(path), + crs=metadata.crs, + shape=metadata.shape, + datetime=metadata.datetime, + ) + return metadata + + +def parse_s1tiling_filename(filename: str) -> dict | None: + """Parse an S1Tiling filename into component fields. + + Returns None if the filename does not match the expected pattern. + """ + m = S1TILING_FILENAME_PATTERN.match(filename) + if not m: + return None + return { + "platform": m.group("platform"), + "tile": m.group("tile"), + "pol": m.group("pol"), + "orbit_dir": m.group("orbit_dir"), + "rel_orbit": m.group("rel_orbit"), + "acq_stamp": m.group("acq_stamp"), + "is_mask": m.group("mask") is not None, + } + + +# ============================================================================= +# Multiscales Layout +# ============================================================================= + + +def compute_multiscales_layout( + native_shape: list[int], + native_transform: list[float], +) -> list[dict]: + """Build the multiscales layout array for all resolution levels.""" + layout: list[dict] = [] + current_shape = native_shape[:] + current_transform = native_transform[:] + + for level_name, parent_name, factor in OVERVIEW_CHAIN: + if parent_name is not None: + current_shape = [ + ceil(current_shape[0] / factor), + ceil(current_shape[1] / factor), + ] + current_transform = [ + current_transform[0] * factor, # a: pixel width + current_transform[1], # b: rotation (0) + current_transform[2], # c: x origin + current_transform[3], # d: rotation (0) + current_transform[4] * factor, # e: pixel height (negative) + current_transform[5], # f: y origin + ] + + entry: dict = { + "asset": level_name, + "spatial:shape": current_shape[:], + "spatial:transform": current_transform[:], + } + if parent_name is None: + entry["transform"] = {"scale": [1.0, 1.0]} + else: + entry["derived_from"] = parent_name + entry["transform"] = { + "scale": [float(factor), float(factor)], + "translation": [0.0, 0.0], + } + + layout.append(entry) + + return layout + + +# ============================================================================= +# Store Creation +# ============================================================================= + + +def _create_spatial_coordinate_arrays( + level_group: zarr.Group, + level_h: int, + level_w: int, + level_transform: list[float], +) -> None: + """Create 1D x and y spatial coordinate arrays at a resolution level.""" + pixel_w = level_transform[0] # a: pixel width + x_origin = level_transform[2] # c: x origin (left edge) + pixel_h = level_transform[4] # e: pixel height (negative) + y_origin = level_transform[5] # f: y origin (top edge) + + x_coords = np.linspace( + x_origin, x_origin + level_w * pixel_w, level_w, endpoint=False, dtype="float64" + ) + y_coords = np.linspace( + y_origin, y_origin + level_h * pixel_h, level_h, endpoint=False, dtype="float64" + ) + + x_arr = level_group.create_array( + "x", + data=x_coords, + chunks=(level_w,), + fill_value=float("nan"), + dimension_names=["x"], + ) + x_arr.attrs.update( + { + "units": "m", + "long_name": "x coordinate of projection", + "standard_name": "projection_x_coordinate", + "_ARRAY_DIMENSIONS": ["x"], + } + ) + + y_arr = level_group.create_array( + "y", + data=y_coords, + chunks=(level_h,), + fill_value=float("nan"), + dimension_names=["y"], + ) + y_arr.attrs.update( + { + "units": "m", + "long_name": "y coordinate of projection", + "standard_name": "projection_y_coordinate", + "_ARRAY_DIMENSIONS": ["y"], + } + ) + + +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. + """ + 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, + "multiscales": { + "layout": layout, + "resampling_method": "average", + }, + "proj:code": metadata.crs, + "spatial:dimensions": ["y", "x"], + "spatial:bbox": metadata.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"] + ) + + # Coordinate variables at native resolution only + r10m = orbit_group["r10m"] + for name, dtype, fill in [ + ("time", "int64", 0), + ("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=" np.ndarray: + """Downsample a 2D array by the given integer factor. + + For average method, handles non-divisible sizes via edge padding. + For nearest method, uses stride-based subsampling. + """ + h, w = data.shape + new_h = ceil(h / factor) + new_w = ceil(w / factor) + + if method == "nearest": + return data[::factor, ::factor][:new_h, :new_w] + + # Average: block mean with edge padding for non-divisible sizes + pad_h = new_h * factor - h + pad_w = new_w * factor - w + if pad_h > 0 or pad_w > 0: + padded = np.pad(data, ((0, pad_h), (0, pad_w)), mode="edge") + else: + padded = data + + reshaped = padded.reshape(new_h, factor, new_w, factor) + if np.issubdtype(data.dtype, np.floating): + return np.nanmean(reshaped, axis=(1, 3)).astype(data.dtype) + return reshaped.mean(axis=(1, 3)).astype(data.dtype) + + +# ============================================================================= +# Acquisition Ingestion +# ============================================================================= + + +def ingest_s1tiling_acquisition( + vv_path: str | Path, + vh_path: str | Path, + border_mask_path: str | Path, + store_path: str | Path, + orbit_direction: str, +) -> int: + """Ingest one S1Tiling acquisition into a GeoZarr V3 store. + + Creates the store if it does not exist, or appends to an existing store. + Returns the time index of the ingested acquisition. + + Parameters + ---------- + vv_path : str or Path + Path to the VV polarisation GeoTIFF. + vh_path : str or Path + Path to the VH polarisation GeoTIFF. + border_mask_path : str or Path + Path to the VV border mask GeoTIFF. + store_path : str or Path + Path to the output Zarr V3 store. + orbit_direction : str + Orbit direction group name (e.g. "ascending", "descending"). + + Returns + ------- + int + The time index (0-based) of the newly ingested acquisition. + + Raises + ------ + FileNotFoundError + If any of the input GeoTIFF paths do not exist. + 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) + store_path = Path(store_path) + + for p in [vv_path, vh_path, border_mask_path]: + if not p.exists(): + raise FileNotFoundError(f"GeoTIFF not found: {p}") + + # Extract metadata from VV file + meta = extract_geotiff_metadata(vv_path) + + log.info( + "Ingesting S1 acquisition", + vv_path=str(vv_path), + orbit_direction=orbit_direction, + ) + + # Create-or-open store + if not store_path.exists(): + root = create_s1_store(store_path, orbit_direction, meta) + 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"] + ) + r10m = orbit_group["r10m"] + for name, dtype, fill in [ + ("time", "int64", 0), + ("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: + """Consolidate metadata at orbit direction and root levels. + + Must be called AFTER all ingestions complete — consolidated metadata + caches array shapes and will become stale if called mid-ingestion. + """ + zarr.consolidate_metadata(str(store_path), path=orbit_direction, zarr_format=3) + zarr.consolidate_metadata(str(store_path), zarr_format=3) + log.info( + "Metadata consolidated", + store_path=str(store_path), + orbit_direction=orbit_direction, + ) + + +# ============================================================================= +# File Discovery +# ============================================================================= + + +def _acq_stamp_from_geotiff(path: str | Path) -> str: + """Resolve a ``YYYYMMDDtHHMMSS`` acquisition stamp from a GeoTIFF's ACQUISITION_DATETIME tag. + + Used when an S1Tiling filename carries a masked multi-frame time (e.g. ``…txxxxxx``), which + has no real timestamp to parse from the name. See #183. + """ + iso = extract_geotiff_metadata(path).datetime # e.g. "2023-01-15T06:12:34" + date_part, time_part = iso.split("T") + return f"{date_part.replace('-', '')}t{time_part.replace(':', '')}" + + +def discover_s1tiling_acquisitions(input_dir: str | Path) -> list[dict]: + """Discover and group S1Tiling GeoTIFF files into acquisition bundles. + + Returns a list of dicts, each with keys: + platform, tile, orbit_dir, rel_orbit, acq_stamp, vv, vh, vv_mask, vh_mask + + Logs warnings for incomplete acquisitions (missing polarisation or mask files). + """ + input_dir = Path(input_dir) + files = sorted(input_dir.glob("*.tif")) + groups: dict[tuple, dict] = {} + + for f in files: + parsed = parse_s1tiling_filename(f.name) + if parsed is None: + continue + + acq_stamp = parsed["acq_stamp"] + if "x" in acq_stamp: + # Multi-frame product: the filename time is masked (…txxxxxx); resolve the real + # stamp from the GeoTIFF ACQUISITION_DATETIME tag so grouping + downstream STAC + # datetime are correct (#183). + acq_stamp = _acq_stamp_from_geotiff(f) + + key = ( + parsed["platform"], + parsed["tile"], + parsed["orbit_dir"], + parsed["rel_orbit"], + acq_stamp, + ) + + if key not in groups: + groups[key] = { + "platform": parsed["platform"], + "tile": parsed["tile"], + "orbit_dir": parsed["orbit_dir"], + "rel_orbit": parsed["rel_orbit"], + "acq_stamp": acq_stamp, + } + + pol = parsed["pol"] + is_mask = parsed["is_mask"] + + if is_mask: + groups[key][f"{pol}_mask"] = f + else: + groups[key][pol] = f + + acquisitions = [] + for key, acq in sorted(groups.items()): + missing = [k for k in ("vv", "vh", "vv_mask", "vh_mask") if k not in acq] + if missing: + log.warning( + "Incomplete acquisition", + key=key, + missing=missing, + ) + acquisitions.append(acq) + + log.info("Discovered acquisitions", count=len(acquisitions), input_dir=str(input_dir)) + return acquisitions + + +# ============================================================================= +# Conditions Ingestion +# ============================================================================= + + +def ingest_s1tiling_conditions( + store_path: str | Path, + orbit_direction: str, + relative_orbit: int, + gamma_area_path: str | Path | None = None, + lia_path: str | Path | None = None, + incidence_angle_path: str | Path | None = None, +) -> None: + """Write time-invariant condition arrays into the conditions group. + + Conditions are per-orbit (not per-acquisition) and have shape (Y, X) only. + The conditions group carries its own proj: and spatial: conventions. + + Parameters + ---------- + store_path : str or Path + Path to an existing Zarr V3 store (must already have the orbit group). + orbit_direction : str + Orbit direction group name (e.g. "ascending", "descending"). + relative_orbit : int + Relative orbit number, used to suffix array names (e.g. 8 → "gamma_area_008"). + gamma_area_path : str, Path, or None + Path to gamma area GeoTIFF. At least one condition path must be provided. + lia_path : str, Path, or None + Path to LIA (sin(LIA)) GeoTIFF. + incidence_angle_path : str, Path, or None + Path to incidence angle GeoTIFF. + + Raises + ------ + ValueError + If no condition paths are provided, or the store/orbit group doesn't exist. + FileNotFoundError + If any provided condition path does not exist. + """ + condition_inputs: list[tuple[str, Path]] = [] + for label, path in [ + ("gamma_area", gamma_area_path), + ("lia", lia_path), + ("incidence_angle", incidence_angle_path), + ]: + if path is not None: + p = Path(path) + if not p.exists(): + raise FileNotFoundError(f"Condition GeoTIFF not found: {p}") + condition_inputs.append((label, p)) + + if not condition_inputs: + raise ValueError("At least one condition path must be provided") + + store_path = Path(store_path) + if not store_path.exists(): + raise ValueError(f"Store does not exist: {store_path}") + + orbit_suffix = f"{relative_orbit:03d}" + + root = zarr.open_group(str(store_path), mode="r+", zarr_format=3) + if orbit_direction not in root: + raise ValueError( + f"Orbit direction '{orbit_direction}' not found in store. " + "Ingest at least one acquisition first." + ) + + orbit = root[orbit_direction] + + # Read reference metadata from the first condition file + ref_label, ref_path = condition_inputs[0] + with 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] + ref_shape = [src.height, src.width] + + # Validate CRS consistency with orbit group + store_crs = dict(orbit.attrs).get("proj:code") + if store_crs and store_crs != ref_crs: + raise ValueError( + f"CRS mismatch: store has {store_crs}, condition GeoTIFF has {ref_crs}" + ) + + # Create or open conditions group + if "conditions" not in orbit: + conditions = orbit.create_group("conditions") + conditions.attrs.update( + { + "proj:code": ref_crs, + "spatial:dimensions": ["y", "x"], + "spatial:transform": ref_transform, + "spatial:shape": ref_shape, + } + ) + log.info("Created conditions group", orbit_direction=orbit_direction) + else: + conditions = orbit["conditions"] + + # Write each condition array + for label, cond_path in condition_inputs: + array_name = f"{label}_{orbit_suffix}" + + with rasterio.open(str(cond_path)) as src: + data = src.read(1).astype(np.float32) + + h, w = data.shape + + if array_name in conditions: + # Overwrite existing array + conditions[array_name][:, :] = data + log.info("Overwrote condition array", array_name=array_name) + else: + arr = conditions.create_array( + array_name, + shape=(h, w), + dtype="float32", + chunks=( + calculate_aligned_chunk_size(h, 512), + calculate_aligned_chunk_size(w, 512), + ), + compressors=zarr.codecs.BloscCodec(cname="zstd", clevel=5), + fill_value=float("nan"), + dimension_names=["y", "x"], + ) + arr[:, :] = data + log.info( + "Wrote condition array", + array_name=array_name, + shape=list(data.shape), + min=float(np.nanmin(data)), + max=float(np.nanmax(data)), + ) + + log.info( + "Conditions ingestion complete", + orbit_direction=orbit_direction, + relative_orbit=orbit_suffix, + arrays=[f"{label}_{orbit_suffix}" for label, _ in condition_inputs], + ) + + +# ============================================================================= +# Conditions File Discovery +# ============================================================================= + + +def discover_s1tiling_conditions(input_dir: str | Path) -> list[dict]: + """Discover S1Tiling condition GeoTIFF files (gamma_area, LIA). + + Returns a list of dicts, each with keys: + tile, orbit, gamma_area (Path), lia (Path or None) + + Groups by (tile, orbit). + """ + input_dir = Path(input_dir) + files = sorted(input_dir.glob("*.tif")) + groups: dict[tuple[str, str], dict] = {} + + for f in files: + m = S1TILING_GAMMA_AREA_PATTERN.match(f.name) + if m: + tile = m.group("tile") + orbit = m.group("orbit") + key = (tile, orbit) + if key not in groups: + groups[key] = {"tile": tile, "orbit": orbit} + groups[key]["gamma_area"] = f + continue + + m = S1TILING_LIA_PATTERN.match(f.name) + if m: + tile = m.group("tile") + orbit = m.group("orbit") + key = (tile, orbit) + if key not in groups: + groups[key] = {"tile": tile, "orbit": orbit} + groups[key]["lia"] = f + + conditions = list(groups.values()) + log.info("Discovered conditions", count=len(conditions), input_dir=str(input_dir)) + return conditions diff --git a/src/eopf_geozarr/data_api/geozarr/common.py b/src/eopf_geozarr/data_api/geozarr/common.py index cbbbd90c..06c4e6c8 100644 --- a/src/eopf_geozarr/data_api/geozarr/common.py +++ b/src/eopf_geozarr/data_api/geozarr/common.py @@ -24,7 +24,7 @@ from pydantic.experimental.missing_sentinel import MISSING from typing_extensions import runtime_checkable -from eopf_geozarr.data_api.geozarr.projjson import ProjJSON # noqa: TC001 +from eopf_geozarr.data_api.geozarr.projjson import ProjJSON if TYPE_CHECKING: from collections.abc import Mapping diff --git a/src/eopf_geozarr/data_api/geozarr/geoproj.py b/src/eopf_geozarr/data_api/geozarr/geoproj.py index c1dd4a7b..d1849aa8 100644 --- a/src/eopf_geozarr/data_api/geozarr/geoproj.py +++ b/src/eopf_geozarr/data_api/geozarr/geoproj.py @@ -8,7 +8,7 @@ from zarr_cm import geo_proj from eopf_geozarr.data_api.geozarr.common import is_none -from eopf_geozarr.data_api.geozarr.projjson import ProjJSON # noqa: TC001 +from eopf_geozarr.data_api.geozarr.projjson import ProjJSON PROJ_UUID = geo_proj.UUID diff --git a/src/eopf_geozarr/data_api/geozarr/multiscales/geozarr.py b/src/eopf_geozarr/data_api/geozarr/multiscales/geozarr.py index e507ae3e..475bd900 100644 --- a/src/eopf_geozarr/data_api/geozarr/multiscales/geozarr.py +++ b/src/eopf_geozarr/data_api/geozarr/multiscales/geozarr.py @@ -8,7 +8,7 @@ # Runtime import (not TYPE_CHECKING): pydantic resolves this annotation when # building MultiscaleGroupAttrs, so the name must exist at runtime. -from zarr_cm import ConventionMetadataObject # noqa: TC002 +from zarr_cm import ConventionMetadataObject from . import tms, zcm diff --git a/src/eopf_geozarr/data_api/geozarr/multiscales/tms.py b/src/eopf_geozarr/data_api/geozarr/multiscales/tms.py index 19525659..484df89c 100644 --- a/src/eopf_geozarr/data_api/geozarr/multiscales/tms.py +++ b/src/eopf_geozarr/data_api/geozarr/multiscales/tms.py @@ -2,7 +2,7 @@ from pydantic import BaseModel -from eopf_geozarr.data_api.geozarr.types import ResamplingMethod # noqa: TC001 +from eopf_geozarr.data_api.geozarr.types import ResamplingMethod class TileMatrix(BaseModel): diff --git a/src/eopf_geozarr/data_api/geozarr/store.py b/src/eopf_geozarr/data_api/geozarr/store.py index 0835ed25..569079dd 100644 --- a/src/eopf_geozarr/data_api/geozarr/store.py +++ b/src/eopf_geozarr/data_api/geozarr/store.py @@ -22,9 +22,7 @@ from eopf_geozarr.data_api.geozarr.multiscales import MultiscaleMeta from eopf_geozarr.data_api.geozarr.multiscales.geozarr import MultiscaleGroupAttrs from eopf_geozarr.data_api.geozarr.multiscales.zcm import ScaleLevel -from eopf_geozarr.data_api.geozarr.projjson import ( - ProjJSON, # noqa: TC001 (runtime use by pydantic) -) +from eopf_geozarr.data_api.geozarr.projjson import ProjJSON class GeoZarrStoreAttrs(BaseModel): diff --git a/src/eopf_geozarr/data_api/s1_rtc.py b/src/eopf_geozarr/data_api/s1_rtc.py new file mode 100644 index 00000000..0654dfc7 --- /dev/null +++ b/src/eopf_geozarr/data_api/s1_rtc.py @@ -0,0 +1,304 @@ +""" +Pydantic-zarr integrated models for Sentinel-1 GRD gamma0T RTC GeoZarr stores. + +Uses the pyz.v3 GroupSpec/ArraySpec with TypedDict members to enforce strict +structure validation — same pattern as s2.py (which uses pyz.v2 for Zarr V2). + +These models validate time-series Zarr V3 stores built from S1Tiling GeoTIFFs +on the Sentinel-2 MGRS grid. This is a *different data product* from the EOPF +L1 GRD models in s1.py — those describe radar-geometry Zarr V2 products. + +Store hierarchy:: + + s1-grd-rtc-{tile}.zarr/ + ├── zarr.json + ├── ascending/ + │ ├── zarr.json # zarr_conventions, multiscales, proj:, spatial: + │ ├── r10m/ # native resolution dataset + │ │ ├── vv/ # (time, Y, X) float32 + │ │ ├── vh/ # (time, Y, X) float32 + │ │ ├── border_mask/ # (time, Y, X) uint8 + │ │ ├── time/ # (time,) int64 datetime + │ │ ├── absolute_orbit/ + │ │ ├── relative_orbit/ + │ │ └── platform/ + │ ├── r20m/ … r720m/ # overview levels (vv, vh, border_mask only) + │ └── conditions/ + │ └── gamma_area_{orbit}/ # (Y, X) float32 + └── descending/ + └── (same structure) +""" + +from __future__ import annotations + +from typing import Any, Literal, Self + +from pydantic import BaseModel, Field, model_validator +from typing_extensions import TypedDict +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.data_api.geozarr.common import DatasetAttrs +from eopf_geozarr.data_api.geozarr.multiscales.zcm import Multiscales +from eopf_geozarr.pyz.v3 import ArraySpec, GroupSpec + +# ============================================================================ +# Constants +# ============================================================================ + +MULTISCALES_UUID = multiscales_cm.UUID +GEO_PROJ_UUID = geo_proj.UUID +SPATIAL_UUID = spatial_cm.UUID + +REQUIRED_CONVENTION_UUIDS = frozenset({MULTISCALES_UUID, GEO_PROJ_UUID, SPATIAL_UUID}) + +ResolutionLevel = Literal["r10m", "r20m", "r60m", "r120m", "r360m", "r720m"] +OrbitDirection = Literal["ascending", "descending"] +Polarisation = Literal["vv", "vh"] + +# ============================================================================ +# Attributes models +# ============================================================================ + + +class S1RtcOrbitGroupAttrs(BaseModel): + """Attributes for an orbit-direction group (ascending or descending). + + Carries the three GeoZarr conventions plus proj:/spatial:/multiscales metadata. + """ + + zarr_conventions: list[dict[str, Any]] + multiscales: Multiscales + proj_code: str = Field(alias="proj:code") + spatial_dimensions: list[str] = Field(alias="spatial:dimensions") + spatial_bbox: list[float] = Field(alias="spatial:bbox") + + model_config = {"extra": "allow", "populate_by_name": True, "serialize_by_alias": True} + + @model_validator(mode="after") + def validate_zarr_conventions(self) -> Self: + """Ensure all three required convention UUIDs are present.""" + present = {c["uuid"] for c in self.zarr_conventions if "uuid" in c} + missing = REQUIRED_CONVENTION_UUIDS - present + if missing: + raise ValueError(f"Missing required zarr_conventions UUIDs: {missing}") + return self + + @model_validator(mode="after") + def validate_spatial_dimensions(self) -> Self: + if self.spatial_dimensions != ["y", "x"]: + raise ValueError( + f"spatial:dimensions must be ['y', 'x'], got {self.spatial_dimensions}" + ) + return self + + @model_validator(mode="after") + def validate_spatial_bbox(self) -> Self: + if len(self.spatial_bbox) != 4: + raise ValueError(f"spatial:bbox must have 4 elements, got {len(self.spatial_bbox)}") + return self + + +class S1RtcResolutionAttrs(BaseModel): + """Attributes for a resolution-level group (r10m, r20m, ...).""" + + spatial_shape: list[int] = Field(alias="spatial:shape") + spatial_transform: list[float] = Field(alias="spatial:transform") + + model_config = {"extra": "allow", "populate_by_name": True, "serialize_by_alias": True} + + @model_validator(mode="after") + def validate_shape(self) -> Self: + if len(self.spatial_shape) != 2: + raise ValueError(f"spatial:shape must have 2 elements, got {len(self.spatial_shape)}") + return self + + @model_validator(mode="after") + def validate_transform(self) -> Self: + if len(self.spatial_transform) != 6: + raise ValueError( + f"spatial:transform must have 6 elements, got {len(self.spatial_transform)}" + ) + return self + + +class S1RtcConditionsAttrs(BaseModel): + """Attributes for the conditions group.""" + + proj_code: str = Field(alias="proj:code") + spatial_dimensions: list[str] = Field(alias="spatial:dimensions") + spatial_transform: list[float] = Field(alias="spatial:transform") + + model_config = {"extra": "allow", "populate_by_name": True, "serialize_by_alias": True} + + +# ============================================================================ +# TypedDict members (same pattern as S2 Sentinel2ResolutionMembers) +# ============================================================================ + + +class S1RtcNativeResolutionMembers(TypedDict, closed=True, total=False): # type: ignore[call-arg] + """Members for the native resolution dataset (r10m). + + Data variables (time, Y, X) plus 1-D coordinate variables (time,). + All fields optional since not all arrays are present during incremental construction. + """ + + vv: ArraySpec[Any] + vh: ArraySpec[Any] + border_mask: ArraySpec[Any] + time: ArraySpec[Any] + absolute_orbit: ArraySpec[Any] + relative_orbit: ArraySpec[Any] + platform: ArraySpec[Any] + + +class S1RtcOverviewResolutionMembers(TypedDict, closed=True, total=False): # type: ignore[call-arg] + """Members for overview resolution datasets (r20m … r720m). + + Only data variables, no coordinate arrays. + """ + + vv: ArraySpec[Any] + vh: ArraySpec[Any] + border_mask: ArraySpec[Any] + + +# ============================================================================ +# Group models (same pattern as S2 Sentinel2ResolutionDataset etc.) +# ============================================================================ + + +class S1RtcNativeResolutionDataset( + GroupSpec[S1RtcResolutionAttrs, S1RtcNativeResolutionMembers] # type: ignore[type-var] +): + """The r10m dataset: data variables + coordinate arrays.""" + + @model_validator(mode="after") + def validate_data_variables(self) -> Self: + """Ensure vv, vh, and border_mask are present.""" + for name in ("vv", "vh", "border_mask"): + if name not in self.members: + raise ValueError(f"Native resolution dataset must contain '{name}' array") + return self + + @property + def vv(self) -> ArraySpec[Any]: + return self.members["vv"] + + @property + def vh(self) -> ArraySpec[Any]: + return self.members["vh"] + + @property + def border_mask(self) -> ArraySpec[Any]: + return self.members["border_mask"] + + +class S1RtcOverviewResolutionDataset( + GroupSpec[S1RtcResolutionAttrs, S1RtcOverviewResolutionMembers] # type: ignore[type-var] +): + """An overview resolution dataset (r20m-r720m): data variables only.""" + + +class S1RtcConditionsGroup(GroupSpec[S1RtcConditionsAttrs, dict[str, ArraySpec[Any]]]): + """Time-invariant condition arrays, keyed by name (e.g. gamma_area_008).""" + + @model_validator(mode="after") + def validate_has_gamma_area(self) -> Self: + """At least one gamma_area_* array should be present.""" + if not any(k.startswith("gamma_area_") for k in self.members): + raise ValueError("Conditions group must contain at least one 'gamma_area_*' array") + return self + + +class S1RtcOrbitGroupMembers(TypedDict, closed=True, total=False): # type: ignore[call-arg] + """Members for an orbit-direction group. + + Contains resolution-level datasets and conditions. + All optional to support incremental store construction. + """ + + r10m: S1RtcNativeResolutionDataset + r20m: S1RtcOverviewResolutionDataset + r60m: S1RtcOverviewResolutionDataset + r120m: S1RtcOverviewResolutionDataset + r360m: S1RtcOverviewResolutionDataset + r720m: S1RtcOverviewResolutionDataset + conditions: S1RtcConditionsGroup + + +class S1RtcOrbitGroup( + GroupSpec[S1RtcOrbitGroupAttrs, S1RtcOrbitGroupMembers] # type: ignore[type-var] +): + """One orbit direction (ascending or descending) with multiscale layout.""" + + @model_validator(mode="after") + def validate_r10m_present(self) -> Self: + if "r10m" not in self.members: + raise ValueError("Orbit group must contain 'r10m' native resolution dataset") + return self + + @property + def r10m(self) -> S1RtcNativeResolutionDataset: + return self.members["r10m"] + + @property + def conditions(self) -> S1RtcConditionsGroup | None: + return self.members.get("conditions") + + def get_resolution(self, level: ResolutionLevel) -> GroupSpec[Any, Any] | None: + """Retrieve a resolution dataset by level name.""" + return self.members.get(level) + + +# ============================================================================ +# Root model (same pattern as S2 Sentinel2Root) +# ============================================================================ + + +class S1RtcRootMembers(TypedDict, closed=True, total=False): # type: ignore[call-arg] + """Members for the root group. At least one orbit direction must be present.""" + + ascending: S1RtcOrbitGroup + descending: S1RtcOrbitGroup + + +class S1RtcRoot(GroupSpec[DatasetAttrs, S1RtcRootMembers]): # type: ignore[type-var] + """Complete S1 GRD RTC GeoZarr V3 hierarchy. + + The hierarchy follows the implementation plan:: + + s1-grd-rtc-{tile}.zarr/ + ├── zarr.json + ├── ascending/ + │ ├── zarr.json # zarr_conventions, multiscales, proj:, spatial: + │ ├── r10m/ + │ │ ├── vv/ # (time, Y, X) float32 + │ │ ├── vh/ # (time, Y, X) float32 + │ │ ├── border_mask/ # (time, Y, X) uint8 + │ │ ├── time/ # (time,) int64 + │ │ ├── absolute_orbit/ + │ │ ├── relative_orbit/ + │ │ └── platform/ + │ ├── r20m/ … r720m/ + │ └── conditions/ + │ └── gamma_area_{orbit}/ + └── descending/ + └── (same) + """ + + @model_validator(mode="after") + def validate_at_least_one_orbit(self) -> Self: + if "ascending" not in self.members and "descending" not in self.members: + raise ValueError("Store must contain at least one orbit group (ascending/descending)") + return self + + @property + def ascending(self) -> S1RtcOrbitGroup | None: + return self.members.get("ascending") + + @property + def descending(self) -> S1RtcOrbitGroup | None: + return self.members.get("descending") diff --git a/tests/_test_data/s1_rtc_examples/s1-grd-rtc-31TCH.json b/tests/_test_data/s1_rtc_examples/s1-grd-rtc-31TCH.json new file mode 100644 index 00000000..de0a4940 --- /dev/null +++ b/tests/_test_data/s1_rtc_examples/s1-grd-rtc-31TCH.json @@ -0,0 +1,1933 @@ +{ + "zarr_format": 3, + "attributes": {}, + "members": { + "descending": { + "zarr_format": 3, + "attributes": { + "zarr_conventions": [ + { + "uuid": "d35379db-88df-4056-af3a-620245f8e347", + "name": "multiscales", + "schema_url": "https://raw.githubusercontent.com/zarr-conventions/multiscales/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-conventions/multiscales/blob/v1/README.md", + "description": "Multiscale layout" + }, + { + "uuid": "f17cb550-5864-4468-aeb7-f3180cfb622f", + "name": "proj:", + "schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md", + "description": "CRS" + }, + { + "uuid": "689b58e2-cf7b-45e0-9fff-9cfc0883d6b4", + "name": "spatial:", + "schema_url": "https://raw.githubusercontent.com/zarr-conventions/spatial/refs/tags/v1/schema.json", + "spec_url": "https://github.com/zarr-conventions/spatial/blob/v1/README.md", + "description": "Spatial coordinates" + } + ], + "multiscales": { + "layout": [ + { + "asset": "r10m", + "spatial:shape": [ + 10980, + 10980 + ], + "spatial:transform": [ + 10.0, + 0.0, + 500000.0, + 0.0, + -10.0, + 5000000.0 + ], + "transform": { + "scale": [ + 1.0, + 1.0 + ] + } + }, + { + "asset": "r20m", + "spatial:shape": [ + 5490, + 5490 + ], + "spatial:transform": [ + 20.0, + 0.0, + 500000.0, + 0.0, + -20.0, + 5000000.0 + ], + "derived_from": "r10m", + "transform": { + "scale": [ + 2.0, + 2.0 + ], + "translation": [ + 0.0, + 0.0 + ] + } + }, + { + "asset": "r60m", + "spatial:shape": [ + 1830, + 1830 + ], + "spatial:transform": [ + 60.0, + 0.0, + 500000.0, + 0.0, + -60.0, + 5000000.0 + ], + "derived_from": "r20m", + "transform": { + "scale": [ + 3.0, + 3.0 + ], + "translation": [ + 0.0, + 0.0 + ] + } + }, + { + "asset": "r120m", + "spatial:shape": [ + 915, + 915 + ], + "spatial:transform": [ + 120.0, + 0.0, + 500000.0, + 0.0, + -120.0, + 5000000.0 + ], + "derived_from": "r60m", + "transform": { + "scale": [ + 2.0, + 2.0 + ], + "translation": [ + 0.0, + 0.0 + ] + } + }, + { + "asset": "r360m", + "spatial:shape": [ + 305, + 305 + ], + "spatial:transform": [ + 360.0, + 0.0, + 500000.0, + 0.0, + -360.0, + 5000000.0 + ], + "derived_from": "r120m", + "transform": { + "scale": [ + 3.0, + 3.0 + ], + "translation": [ + 0.0, + 0.0 + ] + } + }, + { + "asset": "r720m", + "spatial:shape": [ + 153, + 153 + ], + "spatial:transform": [ + 720.0, + 0.0, + 500000.0, + 0.0, + -720.0, + 5000000.0 + ], + "derived_from": "r360m", + "transform": { + "scale": [ + 2.0, + 2.0 + ], + "translation": [ + 0.0, + 0.0 + ] + } + } + ], + "resampling_method": "average" + }, + "proj:code": "EPSG:32631", + "spatial:dimensions": [ + "y", + "x" + ], + "spatial:bbox": [ + 500000.0, + 4890200.0, + 609800.0, + 5000000.0 + ] + }, + "members": { + "r10m": { + "zarr_format": 3, + "attributes": { + "spatial:shape": [ + 10980, + 10980 + ], + "spatial:transform": [ + 10.0, + 0.0, + 500000.0, + 0.0, + -10.0, + 5000000.0 + ] + }, + "members": { + "vv": { + "zarr_format": 3, + "node_type": "array", + "attributes": {}, + "shape": [ + 3, + 10980, + 10980 + ], + "data_type": "float32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 1, + 10980, + 10980 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": "NaN", + "codecs": [ + { + "name": "sharding_indexed", + "configuration": { + "chunk_shape": [ + 1, + 366, + 366 + ], + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "blosc", + "configuration": { + "cname": "zstd", + "clevel": 5 + } + } + ], + "index_codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "crc32c" + } + ], + "index_location": "end" + } + } + ], + "storage_transformers": [], + "dimension_names": [ + "time", + "y", + "x" + ] + }, + "vh": { + "zarr_format": 3, + "node_type": "array", + "attributes": {}, + "shape": [ + 3, + 10980, + 10980 + ], + "data_type": "float32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 1, + 10980, + 10980 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": "NaN", + "codecs": [ + { + "name": "sharding_indexed", + "configuration": { + "chunk_shape": [ + 1, + 366, + 366 + ], + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "blosc", + "configuration": { + "cname": "zstd", + "clevel": 5 + } + } + ], + "index_codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "crc32c" + } + ], + "index_location": "end" + } + } + ], + "storage_transformers": [], + "dimension_names": [ + "time", + "y", + "x" + ] + }, + "border_mask": { + "zarr_format": 3, + "node_type": "array", + "attributes": {}, + "shape": [ + 3, + 10980, + 10980 + ], + "data_type": "uint8", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 1, + 10980, + 10980 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "sharding_indexed", + "configuration": { + "chunk_shape": [ + 1, + 366, + 366 + ], + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "blosc", + "configuration": { + "cname": "zstd", + "clevel": 5 + } + } + ], + "index_codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "crc32c" + } + ], + "index_location": "end" + } + } + ], + "storage_transformers": [], + "dimension_names": [ + "time", + "y", + "x" + ] + }, + "time": { + "zarr_format": 3, + "node_type": "array", + "attributes": {}, + "shape": [ + 3 + ], + "data_type": "int64", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 512 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + } + ], + "storage_transformers": [], + "dimension_names": [ + "time" + ] + }, + "absolute_orbit": { + "zarr_format": 3, + "node_type": "array", + "attributes": {}, + "shape": [ + 3 + ], + "data_type": "int32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 512 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + } + ], + "storage_transformers": [], + "dimension_names": [ + "time" + ] + }, + "relative_orbit": { + "zarr_format": 3, + "node_type": "array", + "attributes": {}, + "shape": [ + 3 + ], + "data_type": "int32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 512 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + } + ], + "storage_transformers": [], + "dimension_names": [ + "time" + ] + }, + "platform": { + "zarr_format": 3, + "node_type": "array", + "attributes": {}, + "shape": [ + 3 + ], + "data_type": " dict[str, object]: @@ -84,6 +85,15 @@ def s2_json_example(request: pytest.FixtureRequest) -> dict[str, object]: return read_json(source_path) +@pytest.fixture(params=s1_rtc_example_json_paths, ids=get_stem) +def s1_rtc_json_example(request: pytest.FixtureRequest) -> dict[str, object]: + """ + A fixture that returns the JSON model of a Sentinel-1 GRD RTC GeoZarr V3 store + """ + source_path: pathlib.Path = request.param + return read_json(source_path) + + @pytest.fixture(params=geozarr_example_paths, ids=get_stem) def s2_geozarr_group_example(request: pytest.FixtureRequest) -> zarr.Group: """ diff --git a/tests/test_data_api/test_s1_rtc.py b/tests/test_data_api/test_s1_rtc.py new file mode 100644 index 00000000..1f531cbc --- /dev/null +++ b/tests/test_data_api/test_s1_rtc.py @@ -0,0 +1,128 @@ +""" +Round-trip and validation tests for Sentinel-1 GRD RTC pydantic-zarr models. + +These tests verify that S1 RTC GeoZarr V3 store metadata can be: +1. Loaded from example JSON data using direct instantiation +2. Validated through Pydantic models +3. Round-tripped without data loss +4. Rejects invalid structures +""" + +from __future__ import annotations + +import copy + +import pytest + +from eopf_geozarr.data_api.s1_rtc import S1RtcRoot + + +def test_s1_rtc_roundtrip(s1_rtc_json_example: dict[str, object]) -> None: + """Test that we can round-trip JSON data without loss.""" + model1 = S1RtcRoot(**s1_rtc_json_example) + dumped = model1.model_dump() + model2 = S1RtcRoot(**dumped) + assert model1.model_dump() == model2.model_dump() + + +def test_s1_rtc_descending_present(s1_rtc_json_example: dict[str, object]) -> None: + """Test that the fixture has a descending orbit group.""" + model = S1RtcRoot(**s1_rtc_json_example) + assert model.descending is not None + assert model.ascending is None + + +def test_s1_rtc_r10m_has_data_arrays(s1_rtc_json_example: dict[str, object]) -> None: + """Test that r10m contains vv, vh, border_mask and coordinate arrays.""" + model = S1RtcRoot(**s1_rtc_json_example) + r10m = model.descending.r10m + assert r10m.vv is not None + assert r10m.vh is not None + assert r10m.border_mask is not None + assert "time" in r10m.members + assert "absolute_orbit" in r10m.members + assert "relative_orbit" in r10m.members + assert "platform" in r10m.members + + +def test_s1_rtc_overview_levels(s1_rtc_json_example: dict[str, object]) -> None: + """Test that overview levels r20m-r720m exist and have vv/vh/border_mask.""" + model = S1RtcRoot(**s1_rtc_json_example) + orbit = model.descending + for level in ("r20m", "r60m", "r120m", "r360m", "r720m"): + group = orbit.get_resolution(level) + assert group is not None, f"Missing overview level {level}" + assert "vv" in group.members + assert "vh" in group.members + assert "border_mask" in group.members + # Overview levels should NOT have coordinate arrays + assert "time" not in group.members + + +def test_s1_rtc_conditions(s1_rtc_json_example: dict[str, object]) -> None: + """Test that conditions group has gamma_area per-orbit arrays.""" + model = S1RtcRoot(**s1_rtc_json_example) + conditions = model.descending.conditions + assert conditions is not None + gamma_keys = [k for k in conditions.members if k.startswith("gamma_area_")] + assert len(gamma_keys) >= 1 + + +def test_s1_rtc_orbit_attrs(s1_rtc_json_example: dict[str, object]) -> None: + """Test that orbit group attributes contain required conventions and metadata.""" + model = S1RtcRoot(**s1_rtc_json_example) + attrs = model.descending.attributes + assert len(attrs.zarr_conventions) == 3 + assert attrs.proj_code.startswith("EPSG:") + assert attrs.spatial_dimensions == ["y", "x"] + assert len(attrs.spatial_bbox) == 4 + assert len(attrs.multiscales.layout) == 6 + assert attrs.multiscales.layout[0].asset == "r10m" + + +def test_s1_rtc_rejects_no_orbit(s1_rtc_json_example: dict[str, object]) -> None: + """Reject a store with no orbit groups.""" + data = copy.deepcopy(s1_rtc_json_example) + data["members"] = {} + with pytest.raises(Exception, match="at least one orbit"): + S1RtcRoot(**data) + + +def test_s1_rtc_rejects_missing_r10m(s1_rtc_json_example: dict[str, object]) -> None: + """Reject an orbit group that lacks r10m.""" + data = copy.deepcopy(s1_rtc_json_example) + del data["members"]["descending"]["members"]["r10m"] + with pytest.raises(Exception, match="r10m"): + S1RtcRoot(**data) + + +def test_s1_rtc_rejects_missing_convention_uuid(s1_rtc_json_example: dict[str, object]) -> None: + """Reject orbit attrs with missing convention UUIDs.""" + data = copy.deepcopy(s1_rtc_json_example) + data["members"]["descending"]["attributes"]["zarr_conventions"] = [ + {"uuid": "fake-uuid", "name": "fake"} + ] + with pytest.raises(Exception, match="Missing required zarr_conventions"): + S1RtcRoot(**data) + + +def test_s1_rtc_rejects_bad_spatial_dimensions(s1_rtc_json_example: dict[str, object]) -> None: + """Reject orbit attrs with wrong spatial:dimensions.""" + data = copy.deepcopy(s1_rtc_json_example) + data["members"]["descending"]["attributes"]["spatial:dimensions"] = ["lat", "lon"] + with pytest.raises(Exception, match="spatial:dimensions"): + S1RtcRoot(**data) + + +def test_s1_rtc_rejects_conditions_without_gamma_area( + s1_rtc_json_example: dict[str, object], +) -> None: + """Reject conditions group with no gamma_area_* arrays.""" + data = copy.deepcopy(s1_rtc_json_example) + cond_members = data["members"]["descending"]["members"]["conditions"]["members"] + # Replace all keys with non-gamma_area names + data["members"]["descending"]["members"]["conditions"]["members"] = { + "some_other": next(iter(cond_members.values())) + } + with pytest.raises(Exception, match="gamma_area"): + S1RtcRoot(**data) diff --git a/tests/test_s1_rtc_ingest.py b/tests/test_s1_rtc_ingest.py new file mode 100644 index 00000000..c5e4b538 --- /dev/null +++ b/tests/test_s1_rtc_ingest.py @@ -0,0 +1,790 @@ +"""Tests for S1 GRD RTC GeoTIFF → GeoZarr V3 ingestion pipeline.""" + +from __future__ import annotations + +from math import ceil +from pathlib import Path + +import numpy as np +import pytest +import rasterio +import xarray as xr +import zarr +from rasterio.transform import from_bounds + +from eopf_geozarr.conversion.s1_ingest import ( + OVERVIEW_CHAIN, + S1TilingMetadata, + _normalise_s1tiling_datetime, + consolidate_s1_store, + create_s1_store, + discover_s1tiling_acquisitions, + discover_s1tiling_conditions, + extract_geotiff_metadata, + ingest_s1tiling_acquisition, + ingest_s1tiling_conditions, + parse_s1tiling_filename, +) + +# ============================================================================= +# Constants +# ============================================================================= + +SIZE = 256 +CRS = "EPSG:32633" +XMIN, YMIN, XMAX, YMAX = 500000.0, 4997440.0, 502560.0, 5000000.0 +TRANSFORM = from_bounds(XMIN, YMIN, XMAX, YMAX, SIZE, SIZE) + +ACQ1_TAGS = { + "ACQUISITION_DATETIME": "2023:01:15T06:12:34Z", + "ORBIT_NUMBER": "47001", + "RELATIVE_ORBIT_NUMBER": "037", + "FLYING_UNIT_CODE": "S1A", + "CALIBRATION": "gamma_naught", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20230115", +} + +ACQ2_TAGS = { + "ACQUISITION_DATETIME": "2023:01:27T06:12:35Z", + "ORBIT_NUMBER": "47177", + "RELATIVE_ORBIT_NUMBER": "037", + "FLYING_UNIT_CODE": "S1A", + "CALIBRATION": "gamma_naught", + "INPUT_S1_IMAGES": "S1A_IW_GRDH_1SDV_20230127", +} + + +# ============================================================================= +# Helpers +# ============================================================================= + + +def _create_synthetic_geotiff( + path: Path, + data: np.ndarray, + crs: str = CRS, + transform: rasterio.transform.Affine | None = None, + tags: dict[str, str] | None = None, +) -> None: + """Write a single-band GeoTIFF with optional metadata tags.""" + if transform is None: + transform = TRANSFORM + with rasterio.open( + str(path), + "w", + driver="GTiff", + height=data.shape[0], + width=data.shape[1], + count=1, + dtype=data.dtype, + crs=crs, + transform=transform, + ) as dst: + if tags: + dst.update_tags(**tags) + dst.write(data, 1) + + +# ============================================================================= +# Fixtures +# ============================================================================= + + +@pytest.fixture +def s1_geotiff_dir(tmp_path: Path) -> Path: + """Create a directory with synthetic S1Tiling GeoTIFFs for 2 acquisitions.""" + rng = np.random.default_rng(42) + + for acq_idx, (stamp, tags) in enumerate( + [("20230115t061234", ACQ1_TAGS), ("20230127t061235", ACQ2_TAGS)] + ): + vv_data = rng.uniform(0.0, 1.0, (SIZE, SIZE)).astype(np.float32) + acq_idx + vh_data = rng.uniform(0.0, 0.5, (SIZE, SIZE)).astype(np.float32) + acq_idx + mask_data = np.ones((SIZE, SIZE), dtype=np.uint8) + mask_data[:10, :] = 0 # border region + + for pol, data in [("vv", vv_data), ("vh", vh_data)]: + fname = f"s1a_32TQM_{pol}_ASC_037_{stamp}_GammaNaughtRTC.tif" + _create_synthetic_geotiff(tmp_path / fname, data, tags=tags) + + mask_fname = f"s1a_32TQM_{pol}_ASC_037_{stamp}_GammaNaughtRTC_BorderMask.tif" + _create_synthetic_geotiff(tmp_path / mask_fname, mask_data, tags=tags) + + return tmp_path + + +@pytest.fixture +def s1_store_path(tmp_path: Path) -> Path: + """Return a clean path for Zarr store output.""" + return tmp_path / "s1-grd-rtc-test.zarr" + + +@pytest.fixture +def single_vv_geotiff(tmp_path: Path) -> Path: + """Create a single VV GeoTIFF with metadata tags.""" + rng = np.random.default_rng(42) + data = rng.uniform(0.0, 1.0, (SIZE, SIZE)).astype(np.float32) + path = tmp_path / "test_vv.tif" + _create_synthetic_geotiff(path, data, tags=ACQ1_TAGS) + return path + + +# ============================================================================= +# Step 9: Metadata extraction tests +# ============================================================================= + + +class TestExtractGeotiffMetadata: + def test_extracts_all_fields(self, single_vv_geotiff: Path) -> None: + meta = extract_geotiff_metadata(single_vv_geotiff) + assert isinstance(meta, S1TilingMetadata) + assert meta.crs == CRS + assert meta.shape == [SIZE, SIZE] + assert len(meta.spatial_transform) == 6 + assert len(meta.bounds) == 4 + assert meta.absolute_orbit == 47001 + assert meta.relative_orbit == 37 + assert meta.platform == "S1A" + assert meta.calibration == "gamma_naught" + + def test_normalises_datetime(self, single_vv_geotiff: Path) -> None: + meta = extract_geotiff_metadata(single_vv_geotiff) + # "2023:01:15T06:12:34Z" → "2023-01-15T06:12:34" + assert meta.datetime == "2023-01-15T06:12:34" + + def test_raises_on_missing_tags(self, tmp_path: Path) -> None: + data = np.zeros((SIZE, SIZE), dtype=np.float32) + path = tmp_path / "no_tags.tif" + _create_synthetic_geotiff(path, data, tags={}) + with pytest.raises(ValueError, match="missing required tags"): + extract_geotiff_metadata(path) + + +class TestNormaliseDatetime: + def test_s1tiling_format(self) -> None: + assert _normalise_s1tiling_datetime("2025:02:10T06:09:20Z") == "2025-02-10T06:09:20" + + def test_already_normalised(self) -> None: + assert _normalise_s1tiling_datetime("2023-01-15T06:12:34") == "2023-01-15T06:12:34" + + +class TestParseFilename: + def test_vv_file(self) -> None: + result = parse_s1tiling_filename("s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC.tif") + assert result is not None + assert result["platform"] == "s1a" + assert result["tile"] == "32TQM" + assert result["pol"] == "vv" + assert result["orbit_dir"] == "ASC" + assert result["rel_orbit"] == "037" + assert result["is_mask"] is False + + def test_mask_file(self) -> None: + result = parse_s1tiling_filename( + "s1a_32TQM_vh_ASC_037_20230115t061234_GammaNaughtRTC_BorderMask.tif" + ) + assert result is not None + assert result["pol"] == "vh" + assert result["is_mask"] is True + + def test_masked_multiframe_time_stamp(self) -> None: + """Multi-frame products carry a masked time (…txxxxxx); the parser must still match so + the file isn't skipped (the real stamp is resolved later from the tag). See #183.""" + result = parse_s1tiling_filename( + "s1a_32TQM_vv_ASC_037_20230115txxxxxx_GammaNaughtRTC.tif" + ) + assert result is not None + assert result["acq_stamp"] == "20230115txxxxxx" + assert result["pol"] == "vv" + + def test_returns_none_for_unknown(self) -> None: + assert parse_s1tiling_filename("random_file.tif") is None + assert parse_s1tiling_filename("not_a_geotiff.txt") is None + + +# ============================================================================= +# Step 10: Store creation tests +# ============================================================================= + + +@pytest.fixture +def sample_metadata(single_vv_geotiff: Path) -> S1TilingMetadata: + """Extract metadata from the single VV fixture.""" + return extract_geotiff_metadata(single_vv_geotiff) + + +class TestCreateStore: + def test_structure(self, s1_store_path: Path, sample_metadata: S1TilingMetadata) -> None: + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + assert "ascending" in root + orbit = root["ascending"] + for level_name, _, _ in OVERVIEW_CHAIN: + assert level_name in orbit, f"Missing level {level_name}" + + def test_conventions(self, s1_store_path: Path, sample_metadata: S1TilingMetadata) -> None: + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + attrs = dict(root["ascending"].attrs) + assert "zarr_conventions" in attrs + conv_names = {c["name"] for c in attrs["zarr_conventions"]} + assert "multiscales" in conv_names + assert "proj:" in conv_names + assert "spatial:" in conv_names + assert attrs["proj:code"] == CRS + assert attrs["spatial:dimensions"] == ["y", "x"] + assert len(attrs["spatial:bbox"]) == 4 + + def test_array_metadata(self, s1_store_path: Path, sample_metadata: S1TilingMetadata) -> None: + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + r10m = root["ascending"]["r10m"] + for arr_name in ["vv", "vh", "border_mask"]: + arr = r10m[arr_name] + assert arr.metadata.dimension_names == ("time", "y", "x") + assert arr.shape[0] == 0 # time axis starts at 0 + assert r10m["vv"].dtype == np.float32 + assert r10m["border_mask"].dtype == np.uint8 + + def test_coordinate_variables( + self, s1_store_path: Path, sample_metadata: S1TilingMetadata + ) -> None: + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + r10m = root["ascending"]["r10m"] + for coord_name in ["time", "absolute_orbit", "relative_orbit", "platform"]: + assert coord_name in r10m, f"Missing coord {coord_name}" + assert r10m[coord_name].shape == (0,) + + def test_overview_shapes(self, s1_store_path: Path, sample_metadata: S1TilingMetadata) -> None: + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + orbit = root["ascending"] + # Verify shape chain follows ceiling division + expected_h, expected_w = SIZE, SIZE + for level_name, _, factor in OVERVIEW_CHAIN: + if factor > 1: + expected_h = ceil(expected_h / factor) + expected_w = ceil(expected_w / factor) + level = orbit[level_name] + arr = level["vv"] + assert arr.shape[1] == expected_h + assert arr.shape[2] == expected_w + + def test_spatial_coordinate_arrays( + self, s1_store_path: Path, sample_metadata: S1TilingMetadata + ) -> None: + """Verify x and y 1D arrays exist at every resolution level.""" + root = create_s1_store(s1_store_path, "ascending", sample_metadata) + orbit = root["ascending"] + for level_name, _, _ in OVERVIEW_CHAIN: + level = orbit[level_name] + for coord in ["x", "y"]: + assert coord in level, f"Missing {coord} at {level_name}" + arr = level[coord] + assert len(arr.shape) == 1 + attrs = dict(arr.attrs) + assert "units" in attrs + assert "standard_name" in attrs + assert "_ARRAY_DIMENSIONS" in attrs + + # Verify x array shape matches level width + level_attrs = dict(level.attrs) + level_h, level_w = level_attrs["spatial:shape"] + assert level["x"].shape[0] == level_w + assert level["y"].shape[0] == level_h + + +# ============================================================================= +# Step 11: Ingestion tests +# ============================================================================= + + +class TestIngestAcquisition: + def _get_acq_paths(self, geotiff_dir: Path, stamp: str) -> tuple[Path, Path, Path]: + """Get VV, VH, border mask paths for a given acquisition stamp.""" + 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_first_acquisition(self, s1_geotiff_dir: Path, s1_store_path: Path) -> None: + vv, vh, mask = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") + idx = ingest_s1tiling_acquisition(vv, vh, mask, s1_store_path, "ascending") + assert idx == 0 + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + assert root["ascending"]["r10m"]["vv"].shape[0] == 1 + + def test_second_acquisition_appends(self, s1_geotiff_dir: Path, s1_store_path: Path) -> None: + 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") + idx = ingest_s1tiling_acquisition(vv2, vh2, mask2, s1_store_path, "ascending") + assert idx == 1 + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + assert root["ascending"]["r10m"]["vv"].shape[0] == 2 + + 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") + + # Read back and compare + with rasterio.open(str(vv)) as src: + expected_vv = src.read(1) + 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) + actual_mask = root["ascending"]["r10m"]["border_mask"][0, :, :] + np.testing.assert_array_equal(actual_mask, expected_mask) + + def test_coordinate_values(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") + + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + r10m = root["ascending"]["r10m"] + assert r10m["absolute_orbit"][0] == 47001 + assert r10m["relative_orbit"][0] == 37 + assert str(r10m["platform"][0]) == "S1A" + + # Verify time is a valid nanosecond timestamp (stored as int64) + time_val = int(r10m["time"][0]) + dt = np.datetime64(time_val, "ns") + assert str(dt).startswith("2023-01-15") + + def test_overview_consistency(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") + + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + orbit = root["ascending"] + expected_h, expected_w = SIZE, SIZE + for level_name, _, factor in OVERVIEW_CHAIN: + if factor > 1: + expected_h = ceil(expected_h / factor) + expected_w = ceil(expected_w / factor) + arr = orbit[level_name]["vv"] + assert arr.shape == (1, expected_h, expected_w), ( + f"Shape mismatch at {level_name}: {arr.shape}" + ) + + def test_rejects_mismatched_crs( + self, s1_geotiff_dir: Path, s1_store_path: Path, tmp_path: Path + ) -> None: + vv1, vh1, mask1 = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(vv1, vh1, mask1, s1_store_path, "ascending") + + # Create a GeoTIFF with different CRS + data = np.ones((SIZE, SIZE), dtype=np.float32) + wrong_crs_dir = tmp_path / "wrong_crs" + wrong_crs_dir.mkdir() + for name, d in [("vv.tif", data), ("vh.tif", data), ("mask.tif", data)]: + _create_synthetic_geotiff(wrong_crs_dir / name, d, crs="EPSG:32632", tags=ACQ1_TAGS) + + with pytest.raises(ValueError, match="CRS mismatch"): + ingest_s1tiling_acquisition( + wrong_crs_dir / "vv.tif", + wrong_crs_dir / "vh.tif", + wrong_crs_dir / "mask.tif", + s1_store_path, + "ascending", + ) + + def test_rejects_mismatched_shape( + self, s1_geotiff_dir: Path, s1_store_path: Path, tmp_path: Path + ) -> None: + vv1, vh1, mask1 = self._get_acq_paths(s1_geotiff_dir, "20230115t061234") + ingest_s1tiling_acquisition(vv1, vh1, mask1, s1_store_path, "ascending") + + # Create GeoTIFFs with different shape + wrong_shape_dir = tmp_path / "wrong_shape" + wrong_shape_dir.mkdir() + small_data = np.ones((128, 128), dtype=np.float32) + small_transform = from_bounds(XMIN, YMIN, XMAX, YMAX, 128, 128) + for name in ["vv.tif", "vh.tif", "mask.tif"]: + _create_synthetic_geotiff( + wrong_shape_dir / name, + small_data, + transform=small_transform, + tags=ACQ1_TAGS, + ) + + with pytest.raises(ValueError, match="Shape mismatch"): + ingest_s1tiling_acquisition( + wrong_shape_dir / "vv.tif", + wrong_shape_dir / "vh.tif", + wrong_shape_dir / "mask.tif", + s1_store_path, + "ascending", + ) + + def test_xarray_roundtrip(self, s1_geotiff_dir: Path, s1_store_path: Path) -> None: + 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, "ascending") + + # Open r10m with xarray + r10m_path = s1_store_path / "ascending" / "r10m" + ds = xr.open_zarr(str(r10m_path)) + assert "vv" in ds + assert ds["vv"].shape[0] == 2 + # Sort by time should work + ds_sorted = ds.sortby("time") + assert ds_sorted["vv"].shape[0] == 2 + + +# ============================================================================= +# Step 12b: Consolidation tests +# ============================================================================= + + +class TestConsolidation: + def test_consolidate_s1_store(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") + consolidate_s1_store(s1_store_path, "ascending") + + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + assert root.metadata.consolidated_metadata is not None + orbit = root["ascending"] + assert orbit.metadata.consolidated_metadata is not None + + def test_consolidate_after_all_ingestions( + self, s1_geotiff_dir: Path, s1_store_path: Path + ) -> None: + 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, "ascending") + consolidate_s1_store(s1_store_path, "ascending") + + # Verify consolidated metadata reflects final shape (2 timesteps) + root = zarr.open_group(str(s1_store_path), mode="r", zarr_format=3) + r10m = root["ascending"]["r10m"] + assert r10m["vv"].shape[0] == 2 + + 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" + mask = geotiff_dir / f"s1a_32TQM_vv_ASC_037_{stamp}_GammaNaughtRTC_BorderMask.tif" + return vv, vh, mask + + +# ============================================================================= +# Step 12: File discovery tests +# ============================================================================= + + +class TestDiscoverAcquisitions: + def test_groups_correctly(self, s1_geotiff_dir: Path) -> None: + acqs = discover_s1tiling_acquisitions(s1_geotiff_dir) + assert len(acqs) == 2 + # Each should have vv, vh, vv_mask, vh_mask + for acq in acqs: + assert "vv" in acq + assert "vh" in acq + assert "vv_mask" in acq + assert "vh_mask" in acq + + def test_warns_on_incomplete(self, tmp_path: Path) -> None: + # Create only VV (no VH, no masks) + data = np.ones((SIZE, SIZE), dtype=np.float32) + fname = "s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC.tif" + _create_synthetic_geotiff(tmp_path / fname, data, tags=ACQ1_TAGS) + + acqs = discover_s1tiling_acquisitions(tmp_path) + assert len(acqs) == 1 + # Should be missing vh, vv_mask, vh_mask + missing = [k for k in ("vh", "vv_mask", "vh_mask") if k not in acqs[0]] + assert len(missing) == 3 + + def test_skips_non_matching(self, tmp_path: Path) -> None: + data = np.ones((SIZE, SIZE), dtype=np.float32) + _create_synthetic_geotiff(tmp_path / "random_file.tif", data, tags=ACQ1_TAGS) + acqs = discover_s1tiling_acquisitions(tmp_path) + assert len(acqs) == 0 + + def test_resolves_masked_multiframe_stamp_from_tag(self, tmp_path: Path) -> None: + """Multi-frame products whose filename time is masked (…txxxxxx) must still be discovered + as a complete acquisition, with acq_stamp resolved from the GeoTIFF ACQUISITION_DATETIME + tag rather than the filename. Regression for #183 (previously returned 0 acquisitions).""" + data = np.ones((SIZE, SIZE), dtype=np.float32) + mask = np.ones((SIZE, SIZE), dtype=np.uint8) + stamp = "20230115txxxxxx" # masked multi-frame time + for pol in ("vv", "vh"): + base = f"s1a_32TQM_{pol}_ASC_037_{stamp}_GammaNaughtRTC" + _create_synthetic_geotiff(tmp_path / f"{base}.tif", data, tags=ACQ1_TAGS) + _create_synthetic_geotiff(tmp_path / f"{base}_BorderMask.tif", mask, tags=ACQ1_TAGS) + + acqs = discover_s1tiling_acquisitions(tmp_path) + + assert len(acqs) == 1 + acq = acqs[0] + # ACQUISITION_DATETIME "2023:01:15T06:12:34Z" -> resolved stamp + assert acq["acq_stamp"] == "20230115t061234" + for k in ("vv", "vh", "vv_mask", "vh_mask"): + assert k in acq + + +# ============================================================================= +# Phase 3: Conditions ingestion tests +# ============================================================================= + + +@pytest.fixture +def s1_store_with_acquisition(s1_geotiff_dir: Path, tmp_path: Path) -> Path: + """Create a Zarr store with one ingested acquisition (prerequisite for conditions).""" + store_path = tmp_path / "s1-grd-rtc-cond.zarr" + vv = s1_geotiff_dir / "s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC.tif" + vh = s1_geotiff_dir / "s1a_32TQM_vh_ASC_037_20230115t061234_GammaNaughtRTC.tif" + mask = s1_geotiff_dir / "s1a_32TQM_vv_ASC_037_20230115t061234_GammaNaughtRTC_BorderMask.tif" + ingest_s1tiling_acquisition(vv, vh, mask, store_path, "ascending") + return store_path + + +@pytest.fixture +def gamma_area_geotiff(tmp_path: Path) -> Path: + """Create a synthetic gamma_area GeoTIFF.""" + rng = np.random.default_rng(99) + data = rng.uniform(0.5, 2.0, (SIZE, SIZE)).astype(np.float32) + path = tmp_path / "GAMMA_AREA_32TQM_037.tif" + _create_synthetic_geotiff(path, data) + return path + + +@pytest.fixture +def lia_geotiff(tmp_path: Path) -> Path: + """Create a synthetic LIA GeoTIFF.""" + rng = np.random.default_rng(100) + data = rng.uniform(0.0, 1.0, (SIZE, SIZE)).astype(np.float32) + path = tmp_path / "sin_LIA_32TQM_037.tif" + _create_synthetic_geotiff(path, data) + return path + + +class TestIngestConditions: + def test_gamma_area_creates_conditions_group( + self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path + ) -> None: + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=gamma_area_geotiff, + ) + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + orbit = root["ascending"] + assert "conditions" in orbit + conditions = orbit["conditions"] + assert "gamma_area_037" in conditions + + def test_conditions_group_attributes( + self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path + ) -> None: + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=gamma_area_geotiff, + ) + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + conditions = root["ascending"]["conditions"] + attrs = dict(conditions.attrs) + assert attrs["proj:code"] == CRS + assert attrs["spatial:dimensions"] == ["y", "x"] + assert len(attrs["spatial:transform"]) == 6 + assert attrs["spatial:shape"] == [SIZE, SIZE] + + def test_gamma_area_array_shape_and_dtype( + self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path + ) -> None: + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=gamma_area_geotiff, + ) + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + arr = root["ascending"]["conditions"]["gamma_area_037"] + assert arr.shape == (SIZE, SIZE) + assert arr.dtype == np.float32 + assert arr.metadata.dimension_names == ("y", "x") + + def test_data_integrity_roundtrip( + self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path + ) -> None: + ingest_s1tiling_conditions( + store_path=s1_store_with_acquisition, + orbit_direction="ascending", + relative_orbit=37, + gamma_area_path=gamma_area_geotiff, + ) + # Read original + with rasterio.open(str(gamma_area_geotiff)) as src: + expected = src.read(1).astype(np.float32) + # Read from Zarr + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + actual = root["ascending"]["conditions"]["gamma_area_037"][:] + np.testing.assert_allclose(actual, expected, rtol=1e-6) + + def test_multiple_conditions( + self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path, lia_geotiff: Path + ) -> None: + 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, + ) + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + conditions = root["ascending"]["conditions"] + assert "gamma_area_037" in conditions + assert "lia_037" in conditions + + def test_multiple_orbits( + self, s1_store_with_acquisition: Path, tmp_path: Path + ) -> None: + """Conditions for different orbits create separate arrays.""" + rng = np.random.default_rng(101) + ga_037 = tmp_path / "GAMMA_AREA_32TQM_037.tif" + ga_110 = tmp_path / "GAMMA_AREA_32TQM_110.tif" + _create_synthetic_geotiff(ga_037, rng.uniform(0.5, 2.0, (SIZE, SIZE)).astype(np.float32)) + _create_synthetic_geotiff(ga_110, rng.uniform(0.5, 2.0, (SIZE, SIZE)).astype(np.float32)) + + ingest_s1tiling_conditions( + s1_store_with_acquisition, "ascending", 37, gamma_area_path=ga_037 + ) + ingest_s1tiling_conditions( + s1_store_with_acquisition, "ascending", 110, gamma_area_path=ga_110 + ) + + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + conditions = root["ascending"]["conditions"] + assert "gamma_area_037" in conditions + assert "gamma_area_110" in conditions + + def test_overwrite_existing_condition( + self, s1_store_with_acquisition: Path, tmp_path: Path + ) -> None: + """Writing the same condition array twice overwrites data.""" + rng = np.random.default_rng(102) + ga_path = tmp_path / "GAMMA_AREA_32TQM_037.tif" + + data_v1 = np.ones((SIZE, SIZE), dtype=np.float32) + _create_synthetic_geotiff(ga_path, data_v1) + ingest_s1tiling_conditions( + s1_store_with_acquisition, "ascending", 37, gamma_area_path=ga_path + ) + + data_v2 = np.full((SIZE, SIZE), 2.0, dtype=np.float32) + _create_synthetic_geotiff(ga_path, data_v2) + ingest_s1tiling_conditions( + s1_store_with_acquisition, "ascending", 37, gamma_area_path=ga_path + ) + + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + actual = root["ascending"]["conditions"]["gamma_area_037"][:] + np.testing.assert_allclose(actual, data_v2, rtol=1e-6) + + def test_raises_no_conditions_provided( + self, s1_store_with_acquisition: Path + ) -> None: + with pytest.raises(ValueError, match="At least one condition"): + ingest_s1tiling_conditions( + s1_store_with_acquisition, "ascending", 37 + ) + + def test_raises_store_not_exists(self, tmp_path: Path, gamma_area_geotiff: Path) -> None: + with pytest.raises(ValueError, match="Store does not exist"): + ingest_s1tiling_conditions( + tmp_path / "nonexistent.zarr", "ascending", 37, + gamma_area_path=gamma_area_geotiff, + ) + + def test_raises_orbit_not_exists( + self, tmp_path: Path, gamma_area_geotiff: Path + ) -> None: + """Raise if the orbit group hasn't been created yet.""" + # Create minimal empty store + store_path = tmp_path / "empty-store.zarr" + zarr.open_group(str(store_path), mode="w-", zarr_format=3) + with pytest.raises(ValueError, match="not found in store"): + ingest_s1tiling_conditions( + store_path, "ascending", 37, gamma_area_path=gamma_area_geotiff + ) + + def test_raises_file_not_found(self, s1_store_with_acquisition: Path) -> None: + with pytest.raises(FileNotFoundError): + ingest_s1tiling_conditions( + s1_store_with_acquisition, "ascending", 37, + gamma_area_path="/nonexistent/gamma_area.tif", + ) + + def test_consolidation_includes_conditions( + self, s1_store_with_acquisition: Path, gamma_area_geotiff: Path + ) -> None: + """Consolidation after conditions ingestion includes the conditions group.""" + ingest_s1tiling_conditions( + s1_store_with_acquisition, "ascending", 37, gamma_area_path=gamma_area_geotiff + ) + consolidate_s1_store(s1_store_with_acquisition, "ascending") + + root = zarr.open_group(str(s1_store_with_acquisition), mode="r", zarr_format=3) + assert root.metadata.consolidated_metadata is not None + orbit = root["ascending"] + assert orbit.metadata.consolidated_metadata is not None + # Conditions group should be accessible through consolidated metadata + assert "conditions" in orbit + assert "gamma_area_037" in orbit["conditions"] + + +# ============================================================================= +# Phase 3: Conditions file discovery tests +# ============================================================================= + + +class TestDiscoverConditions: + def test_discovers_gamma_area(self, tmp_path: Path) -> None: + data = np.ones((SIZE, SIZE), dtype=np.float32) + _create_synthetic_geotiff(tmp_path / "GAMMA_AREA_32TQM_037.tif", data) + _create_synthetic_geotiff(tmp_path / "GAMMA_AREA_32TQM_110.tif", data) + + conditions = discover_s1tiling_conditions(tmp_path) + assert len(conditions) == 2 + orbits = {c["orbit"] for c in conditions} + assert orbits == {"037", "110"} + for c in conditions: + assert "gamma_area" in c + assert c["tile"] == "32TQM" + + def test_discovers_lia(self, tmp_path: Path) -> None: + data = np.ones((SIZE, SIZE), dtype=np.float32) + _create_synthetic_geotiff(tmp_path / "sin_LIA_32TQM_037.tif", data) + + conditions = discover_s1tiling_conditions(tmp_path) + assert len(conditions) == 1 + assert "lia" in conditions[0] + + def test_groups_gamma_area_and_lia(self, tmp_path: Path) -> None: + """Gamma area and LIA for the same tile/orbit are grouped together.""" + data = np.ones((SIZE, SIZE), dtype=np.float32) + _create_synthetic_geotiff(tmp_path / "GAMMA_AREA_32TQM_037.tif", data) + _create_synthetic_geotiff(tmp_path / "sin_LIA_32TQM_037.tif", data) + + conditions = discover_s1tiling_conditions(tmp_path) + assert len(conditions) == 1 + assert "gamma_area" in conditions[0] + assert "lia" in conditions[0] + assert conditions[0]["tile"] == "32TQM" + assert conditions[0]["orbit"] == "037" + + def test_skips_non_matching(self, tmp_path: Path) -> None: + data = np.ones((SIZE, SIZE), dtype=np.float32) + _create_synthetic_geotiff(tmp_path / "random_file.tif", data) + conditions = discover_s1tiling_conditions(tmp_path) + assert len(conditions) == 0 + + def test_empty_directory(self, tmp_path: Path) -> None: + conditions = discover_s1tiling_conditions(tmp_path) + assert len(conditions) == 0