From 3ba9e853564649f03208bfa8f620aaedae798e1c Mon Sep 17 00:00:00 2001 From: ilias <121384694+iliasreg@users.noreply.github.com> Date: Thu, 5 Nov 2026 14:14:29 +0100 Subject: [PATCH 1/3] Read the SAR metadata from the CDSE daily parquet files A recipe with sar_parquet_root takes the imagette footprints, times and centres from the daily parquet files (WVImagettes plus the footprint polygons) and opens no SAFE. Dateline halves are merged in either order, and polygon k is paired with imagette k. SAFEs without an imagette list, SAFEs whose polygons and imagettes disagree, and imagettes whose file is not in the archive are named in the manifest discovery and logged at each run. Adds pyarrow as a dependency. --- environment.yml | 1 + pyproject.toml | 2 +- src/soba_matchups/config.py | 14 +- src/soba_matchups/pipeline.py | 55 +++++-- src/soba_matchups/sar_parquet.py | 268 +++++++++++++++++++++++++++++++ tests/data/parquet_rows.json | 30 ++++ tests/test_sar_parquet.py | 228 ++++++++++++++++++++++++++ 7 files changed, 584 insertions(+), 14 deletions(-) create mode 100644 src/soba_matchups/sar_parquet.py create mode 100644 tests/data/parquet_rows.json create mode 100644 tests/test_sar_parquet.py diff --git a/environment.yml b/environment.yml index 615ab61..24973eb 100644 --- a/environment.yml +++ b/environment.yml @@ -15,3 +15,4 @@ dependencies: - xarray=2026.9.0 - shapely=2.2.0 - netcdf4=1.7.4 + - pyarrow=26.0.0 diff --git a/pyproject.toml b/pyproject.toml index eacc773..0771302 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,7 +11,7 @@ license = "MIT" license-files = ["LICENSE"] requires-python = ">=3.11" dependencies = ["numpy>=1.26,<3", "xarray>=2024.1", "shapely>=2", - "netCDF4>=1.7,<2"] + "netCDF4>=1.7,<2", "pyarrow>=15"] [project.scripts] soba_matchups = "soba_matchups.cli:main" diff --git a/src/soba_matchups/config.py b/src/soba_matchups/config.py index 3e3e6ae..01a56f3 100644 --- a/src/soba_matchups/config.py +++ b/src/soba_matchups/config.py @@ -43,6 +43,7 @@ class Recipe: instrument: str platform: str resolution: str + sar_parquet_root: Path | None = None def resolved(self): """Return the recipe as plain values that can be written to JSON.""" @@ -57,8 +58,8 @@ def load_recipe(path): with path.open("rb") as stream: data = tomllib.load(stream) # Reject unknown keys and missing required keys. - required = set(Recipe.__dataclass_fields__) - {"path"} - unknown = set(data) - required + required = set(Recipe.__dataclass_fields__) - {"path", "sar_parquet_root"} + unknown = set(data) - required - {"sar_parquet_root"} if unknown or required - set(data): raise ValueError( f"Unknown recipe keys: {sorted(unknown)}; missing keys: {sorted(required - set(data))}" @@ -99,12 +100,17 @@ def load_recipe(path): ): raise ValueError("sar_missions must contain distinct supported Sentinel-1 missions") # Make every path absolute, relative to the recipe file. - for key in ["reference_root", "sar_archive_root", "output_dir"]: + paths = ["reference_root", "sar_archive_root", "output_dir"] + inputs = ["reference_root", "sar_archive_root"] + if "sar_parquet_root" in data: + paths.append("sar_parquet_root") + inputs.append("sar_parquet_root") + for key in paths: if not isinstance(data[key], str) or not data[key].strip(): raise ValueError(f"{key} must be a nonempty path") data[key] = (path.parent / data[key]).resolve() # Inputs must exist, and results must never be written inside them. - for key in ["reference_root", "sar_archive_root"]: + for key in inputs: if not data[key].is_dir(): raise ValueError(f"{key} directory does not exist: {data[key]}") if data["output_dir"].is_relative_to(data[key]): diff --git a/src/soba_matchups/pipeline.py b/src/soba_matchups/pipeline.py index 7e632c6..94b3482 100644 --- a/src/soba_matchups/pipeline.py +++ b/src/soba_matchups/pipeline.py @@ -10,6 +10,12 @@ from soba_matchups.adapters import alti, scat, swot from soba_matchups.config import REFERENCE_PRODUCTS from soba_matchups.output import file_identity, imagette_item, publish +from soba_matchups.sar_parquet import ( + IMAGETTE_MAPPING, + TIME_SOURCE, + discover_parquet, + verify_parquet_inputs, +) from soba_matchups.sar_metadata import ( OCN_HEADER_PRECISION_SECONDS, decode_safe_name, @@ -18,6 +24,22 @@ ) +# What each SAR source can and cannot do, copied into the manifest of every run. +ARCHIVE_LIMITATIONS = [ + "Requires extracted L2 OCN SAFEs with manifest.safe and every listed measurement file.", + "Missing/ambiguous native frame associations fail without suffix fallback.", + "Reads OSW geolocation and UTC attributes, not wind/spectral arrays; no network.", + "Archive inventory only; missing mission/day directories are listed in discovery.", + "SAR file identities use size/mtime; manifests and references use SHA256.", +] +PARQUET_LIMITATIONS = [ + "SAR metadata comes from the daily parquet files; no SAFE is opened, nothing is fetched.", + "Imagette centres are polygon centroids; polygon k is paired with imagette number k.", + "SAFEs with no imagette list, or with polygons and imagettes that disagree, are left out.", + "Imagettes whose file is not in the archive are matched and listed as unavailable.", + "Parquet files and references use SHA256.", +] + # One adapter per instrument: it finds the files, reads them and names the point fields. ADAPTERS = {"ascat": scat, "hscat": scat, "karin": swot, "alti": alti} @@ -75,7 +97,17 @@ def iso(t): # Search window: reference coverage widened by the time tolerance and header rounding. start = iso(min(ref["coverage_start"] for ref in refs) - delta - precision) end = iso(max(ref["coverage_end"] for ref in refs) + delta + precision) - products, discovery = discover_local(start, end, recipe.sar_missions, recipe.sar_archive_root) + from_parquet = recipe.sar_parquet_root is not None + if from_parquet: + products, discovery = discover_parquet( + start, end, recipe.sar_missions, recipe.sar_parquet_root, recipe.sar_archive_root + ) + verify_inputs = verify_parquet_inputs + else: + products, discovery = discover_local( + start, end, recipe.sar_missions, recipe.sar_archive_root + ) + verify_inputs = verify_local_inputs items, diagnostics = [], [] seen_points = set() imagettes_processed = duplicate_points = safes_without_candidates = 0 @@ -189,7 +221,7 @@ def number(x, integer=False): imagette_item(product, imagette, fields, groups, recipe.time_tolerance_minutes) ) # Confirm no input changed while matching, before anything is published. - verify_local_inputs(discovery) + verify_inputs(discovery) for ref in refs: if file_identity(ref["path"]) != ref["identity"]: raise ValueError(f"Input changed during processing: {ref['path']}") @@ -227,6 +259,9 @@ def number(x, integer=False): "safe_candidate_cell_rows": candidate_cells, "imagettes_processed": imagettes_processed, "imagettes_skipped": sum(len(p["SkippedImagettes"]) for p in products), + "safes_without_imagettes": len(discovery.get("safes_without_imagettes", [])), + "safes_skipped": len(discovery.get("safes_skipped", [])), + "imagettes_unavailable": len(discovery.get("imagettes_unavailable", [])), "stac_items": len(items), "matchup_imagette_cell_rows": len(seen_points), "duplicate_associations_discarded": duplicate_points, @@ -235,12 +270,18 @@ def number(x, integer=False): "source_code": [file_identity(p) for p in sorted(Path(__file__).parent.rglob("*.py"))], "rules": { "policy": adapter.POLICY, - "imagette_mapping": "local_ocn_geolocation_to_manifest_frame", + "imagette_mapping": ( + IMAGETTE_MAPPING if from_parquet else "local_ocn_geolocation_to_manifest_frame" + ), "space": getattr( adapter, "SPACE", "strict individual imagette footprint point containment" ), "time": "strict imagette start - tolerance < cell time < imagette end + tolerance", - "time_source": "OCN firstMeasurementTime / lastMeasurementTime UTC attributes", + "time_source": ( + TIME_SOURCE + if from_parquet + else "OCN firstMeasurementTime / lastMeasurementTime UTC attributes" + ), "candidate_padding_seconds": OCN_HEADER_PRECISION_SECONDS, "header_precision": "whole-second NC headers use outward-rounded manifest bounds", "time_tolerance_minutes": recipe.time_tolerance_minutes, @@ -252,11 +293,7 @@ def number(x, integer=False): "ref_lat": "degrees_north", }, "limitations": [ - "Requires extracted L2 OCN SAFEs with manifest.safe and all listed measurement files.", - "Missing/ambiguous native frame associations fail without suffix fallback.", - "Reads OSW geolocation and UTC attributes, not wind/spectral arrays; no network.", - "Archive inventory only; missing mission/day directories are listed in discovery.", - "SAR file identities use size/mtime; manifests and references use SHA256.", + *(PARQUET_LIMITATIONS if from_parquet else ARCHIVE_LIMITATIONS), "sarmatchup:* properties are project-specific, not a registered STAC extension.", "No sea-percent, polar, quality or ice filter.", "Tolerance is explicitly configured, not scientifically calibrated.", diff --git a/src/soba_matchups/sar_parquet.py b/src/soba_matchups/sar_parquet.py new file mode 100644 index 0000000..4b7753f --- /dev/null +++ b/src/soba_matchups/sar_parquet.py @@ -0,0 +1,268 @@ +"""SAR metadata from the CDSE daily parquet files, with no SAFE opened and nothing fetched. + +A daily file holds one row per Sentinel-1 product. For a WV OCN product the footprint is one +polygon per imagette, and the column WVImagettes lists the imagette files (name, number, times). +Polygon k belongs to the imagette numbered k, once the two halves of an imagette cut at the +dateline are merged. The archive is asked only which imagette files exist, by name. +""" + +from datetime import datetime, timedelta +import hashlib +import json +import logging +import math +import re +from pathlib import Path + +import numpy as np +import pyarrow.parquet as pq +from shapely.geometry import Polygon, mapping, shape +from shapely.ops import unary_union + +from .output import file_identity +from .sar_metadata import _utc + +logger = logging.getLogger(__name__) +IMAGETTE_MAPPING = "parquet_polygon_number_equals_file_suffix_after_dateline_merge" +TIME_SOURCE = "WVImagettes ContentDate (whole seconds from the file name)" +COLUMNS = ["Name", "Id", "ContentDate", "GeoFootprint", "WVImagettes"] +# A vertex this close to +-180 degrees lies on the dateline cut. +CUT_EDGE_TOLERANCE = 1e-6 +# Largest latitude difference, in degrees, between the two ends of a cut shared by two halves. +CUT_LATITUDE_TOLERANCE = 1e-3 +_WV_OCN = re.compile(r"S1[ABCD]_WV_OCN__2[A-Z]") + + +def _rings(footprint): + """Return the polygons of a GeoJSON footprint as lists of (lon, lat), closing vertex dropped.""" + geometry = json.loads(footprint) + if geometry["type"] == "Polygon": + polygons = [geometry["coordinates"]] + elif geometry["type"] == "MultiPolygon": + polygons = geometry["coordinates"] + else: + raise ValueError(f"Unexpected footprint type {geometry['type']}") + return [[tuple(vertex) for vertex in polygon[0][:-1]] for polygon in polygons] + + +def _cut_latitudes(ring, side): + """Return the sorted latitudes of the vertices on the dateline cut (side is 180 or -180).""" + return sorted(y for x, y in ring if abs(x - side) < CUT_EDGE_TOLERANCE) + + +def _are_halves(first, second): + """Tell whether two polygons are the east and west halves of one imagette, in either order.""" + for east, west in ((first, second), (second, first)): + a, b = _cut_latitudes(east, 180), _cut_latitudes(west, -180) + if a and b and abs(a[0] - b[0]) < CUT_LATITUDE_TOLERANCE: + if abs(a[-1] - b[-1]) < CUT_LATITUDE_TOLERANCE: + return True + return False + + +def merge_dateline(rings): + """Group the polygons that are halves of one imagette; return lists of polygon positions. + + Two neighbours in the footprint are merged when one has vertices on the +180 cut, the other + on the -180 cut, and the latitudes at the two ends of the cut agree. Either may come first. + """ + groups, position = [], 0 + while position < len(rings): + if position + 1 < len(rings) and _are_halves(rings[position], rings[position + 1]): + groups.append([position, position + 1]) + position += 2 + else: + groups.append([position]) + position += 1 + return groups + + +def centre(rings): + """Return the area centroid [lon, lat] of an imagette, given its polygon or its two halves.""" + lat0 = float(np.mean([y for ring in rings for _, y in ring])) + lon0 = rings[0][0][0] + scale = math.cos(math.radians(lat0)) + # A local plane in degrees keeps a dateline pair contiguous; longitudes unwrap around lon0. + planar = [ + Polygon([((((x - lon0 + 180) % 360) - 180) * scale, y - lat0) for x, y in ring]) + for ring in rings + ] + union = unary_union(planar) + if union.is_empty or union.area == 0: + raise ValueError("Empty imagette footprint") + middle = union.centroid + return [float((lon0 + middle.x / scale + 180) % 360 - 180), float(lat0 + middle.y)] + + +def _geometry(rings): + """Return the GeoJSON geometry of an imagette: one Polygon, or a MultiPolygon for two halves.""" + shapes = [Polygon(ring) for ring in rings] + if any(not s.is_valid or s.is_empty for s in shapes): + raise ValueError("Invalid imagette polygon") + return mapping(shapes[0] if len(shapes) == 1 else unary_union(shapes)) + + +def _safe_folder(archive_root, name): + """Return the archive folder of a WV OCN SAFE (sentinel-1x/L2/WV/_WV_OCN__2S/year/day).""" + mission = name[:3] + start = datetime.strptime(name.split("_")[5], "%Y%m%dT%H%M%S") + return ( + archive_root + / ("sentinel-1" + mission[-1].lower()) + / "L2" + / "WV" + / (mission + "_WV_OCN__2S") + / start.strftime("%Y/%j") + / name + ) + + +def _read_day(path): + """Return the rows of one daily parquet file as dicts, reading only the columns used.""" + missing = [c for c in COLUMNS if c not in pq.read_schema(path).names] + if missing: + raise ValueError(f"Daily parquet file lacks the columns {missing}: {path}") + return pq.read_table(path, columns=COLUMNS).to_pylist() + + +def _product(row, archive_root): + """Build the product record of one parquet row, with its imagettes and their availability.""" + name = row["Name"] + listed = json.loads(row["WVImagettes"]) + content = json.loads(row["ContentDate"]) + _utc(content["Start"]), _utc(content["End"]) + rings = _rings(row["GeoFootprint"]) + groups = merge_dateline(rings) + numbers = sorted(item["index"] for item in listed) + if len(groups) != len(listed): + raise ValueError( + f"{len(groups)} footprint polygons after the dateline merge, {len(listed)} imagettes" + ) + if numbers != list(range(1, len(numbers) + 1)): + raise ValueError("Imagette numbers are not 1..N") + product_id = "local_" + hashlib.sha256(name.encode()).hexdigest() + folder = _safe_folder(archive_root, name) + present = ( + {p.name for p in (folder / "measurement").iterdir()} + if (folder / "measurement").is_dir() + else None + ) + imagettes, unavailable = [], [] + for item in sorted(listed, key=lambda item: item["index"]): + group = groups[item["index"] - 1] + member = [rings[i] for i in group] + # The same identifier as the local archive source gives, so both sources agree on Item ids. + relative = "measurement/" + item["Name"] + imagettes.append( + { + "Id": "local_" + hashlib.sha256((product_id + "/" + relative).encode()).hexdigest(), + "Name": item["Name"], + "index": item["index"], + "ContentDate": item["ContentDate"], + "GeoFootprint": _geometry(member), + "Centre": centre(member), + "href": (folder / relative).as_uri(), + "available": present is not None and item["Name"] in present, + } + ) + if not imagettes[-1]["available"]: + reason = "SAFE not in the archive" if present is None else "file not in the archive" + unavailable.append( + {"safe": name, "index": item["index"], "name": item["Name"], "reason": reason} + ) + return { + "Id": product_id, + "Name": name, + "ContentDate": content, + "GeoFootprint": mapping(unary_union([shape(i["GeoFootprint"]) for i in imagettes])), + "Imagettes": imagettes, + "SkippedImagettes": [], + }, unavailable + + +def discover_parquet( + start: str, end: str, missions: tuple[str, ...], parquet_root: Path, sar_archive_root: Path +) -> tuple[list[dict], dict]: + """Read the daily parquet files of the covered days and return the WV OCN products to match. + + Products with no imagette list, and products whose polygons and imagettes disagree, are left + out and named in the returned provenance. Imagette files absent from the archive are kept, + flagged unavailable and named there too. Coverage must already include the time tolerance. + """ + first, last = _utc(start), _utc(end) + if first >= last: + raise ValueError("start must precede end") + provenance = { + "source": "cdse_daily_parquet", + "parquet_inputs": [], + "scanned_files": [], + "missing_files": [], + "coverage": {"start": start, "end": end}, + "missions": sorted(set(missions)), + "sar_parquet_root": str(parquet_root), + "sar_archive_root": str(sar_archive_root), + "coverage_scope": "products listed in the selected daily parquet files", + "safes_without_imagettes": [], + "safes_skipped": [], + "imagettes_unavailable": [], + } + products, seen, count = [], set(), 0 + day = first.date() - timedelta(days=1) + while day <= last.date(): + path = parquet_root / f"s1_{day.isoformat()}.parquet" + provenance["scanned_files"].append(str(path)) + if not path.is_file(): + provenance["missing_files"].append(str(path)) + else: + provenance["parquet_inputs"].append(file_identity(path)) + for row in _read_day(path): + name = row["Name"] + if not _WV_OCN.match(name) or name[:3] not in missions or name in seen: + continue + seen.add(name) + count += 1 + content = json.loads(row["ContentDate"]) + if not (_utc(content["Start"]) < last and _utc(content["End"]) > first): + continue + if row["WVImagettes"] is None: + provenance["safes_without_imagettes"].append( + {"safe": name, "reason": "no imagette list in the parquet file"} + ) + continue + try: + product, unavailable = _product(row, sar_archive_root) + except (ValueError, KeyError, TypeError) as exc: + provenance["safes_skipped"].append({"safe": name, "reason": str(exc)}) + continue + provenance["imagettes_unavailable"].extend(unavailable) + products.append(product) + day += timedelta(days=1) + if not provenance["parquet_inputs"]: + raise ValueError(f"No daily parquet files found under: {parquet_root}") + provenance["safes_found"] = count + provenance["products_selected"] = len(products) + # Everything left out or unusable is named in the manifest and summarised here at each run. + for key, text in [ + ("safes_without_imagettes", "SAFEs without an imagette list"), + ("safes_skipped", "SAFEs skipped (polygons and imagettes disagree)"), + ("imagettes_unavailable", "imagettes whose file is not in the archive"), + ("missing_files", "daily parquet files missing"), + ]: + if provenance[key]: + logger.warning( + "%d %s; listed in the manifest discovery.%s", len(provenance[key]), text, key + ) + return sorted(products, key=lambda p: (p["ContentDate"]["Start"], p["Name"])), provenance + + +def verify_parquet_inputs(provenance: dict) -> None: + """Fail closed when any daily parquet file changed during the run.""" + if provenance.get("source") != "cdse_daily_parquet": + raise ValueError("Expected cdse_daily_parquet provenance") + for identity in provenance["parquet_inputs"]: + try: + unchanged = file_identity(Path(identity["path"])) == identity + except OSError as exc: + raise ValueError(f"Parquet input changed or missing: {identity['path']}") from exc + if not unchanged: + raise ValueError(f"Parquet input changed: {identity['path']}") diff --git a/tests/data/parquet_rows.json b/tests/data/parquet_rows.json new file mode 100644 index 0000000..22bca94 --- /dev/null +++ b/tests/data/parquet_rows.json @@ -0,0 +1,30 @@ +[ + { + "Name": "S1A_WV_OCN__2SSV_20241221T173322_20241221T173508_057093_0704A8_A00B.SAFE", + "Id": "f5f5f50d-c929-412c-85d1-800b88fd2caa", + "ContentDate": "{\"Start\": \"2024-12-21T17:33:22.925957Z\", \"End\": \"2024-12-21T17:35:08.386423Z\"}", + "GeoFootprint": "{\"type\": \"MultiPolygon\", \"coordinates\": [[[[-178.566788, -19.907743], [-178.523682, -19.727993], [-178.719955, -19.684353], [-178.763367, -19.864035], [-178.566788, -19.907743]]], [[[179.483902, -20.393089], [179.528595, -20.219107], [179.337341, -20.174488], [179.292419, -20.348402], [179.483902, -20.393089]]], [[[-179.003754, -21.670485], [-178.959763, -21.490974], [-179.15863, -21.446774], [-179.202957, -21.626211], [-179.003754, -21.670485]]], [[[179.016479, -22.149361], [179.062347, -21.97551], [178.868927, -21.930225], [178.8228, -22.104004], [179.016479, -22.149361]]], [[[-179.450989, -23.432484], [-179.405914, -23.252865], [-179.607635, -23.208067], [-179.653076, -23.387611], [-179.450989, -23.432484]]], [[[178.540848, -23.904806], [178.587936, -23.731356], [178.391983, -23.685343], [178.344604, -23.858717], [178.540848, -23.904806]]], [[[-179.90918, -25.193073], [-179.862991, -25.013683], [-180.0, -24.983298454582286], [-180.0, -25.17293615461778], [-179.90918, -25.193073]]], [[[180.0, -25.17293615461778], [180.0, -24.983298454582286], [179.932144, -24.96825], [179.885559, -25.147562], [180.0, -25.17293615461778]]], [[[178.0513, -25.658716], [178.099838, -25.485281], [177.901108, -25.438499], [177.852249, -25.611855], [178.0513, -25.658716]]]]}", + "WVImagettes": "[{\"Name\": \"s1a-wv1-ocn-vv-20241221t173322-20241221t173325-057093-0704a8-001.nc\", \"index\": 1, \"ContentDate\": {\"Start\": \"2024-12-21T17:33:22Z\", \"End\": \"2024-12-21T17:33:25Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t173337-20241221t173340-057093-0704a8-002.nc\", \"index\": 2, \"ContentDate\": {\"Start\": \"2024-12-21T17:33:37Z\", \"End\": \"2024-12-21T17:33:40Z\"}}, {\"Name\": \"s1a-wv1-ocn-vv-20241221t173352-20241221t173355-057093-0704a8-003.nc\", \"index\": 3, \"ContentDate\": {\"Start\": \"2024-12-21T17:33:52Z\", \"End\": \"2024-12-21T17:33:55Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t173406-20241221t173409-057093-0704a8-004.nc\", \"index\": 4, \"ContentDate\": {\"Start\": \"2024-12-21T17:34:06Z\", \"End\": \"2024-12-21T17:34:09Z\"}}, {\"Name\": \"s1a-wv1-ocn-vv-20241221t173421-20241221t173424-057093-0704a8-005.nc\", \"index\": 5, \"ContentDate\": {\"Start\": \"2024-12-21T17:34:21Z\", \"End\": \"2024-12-21T17:34:24Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t173436-20241221t173439-057093-0704a8-006.nc\", \"index\": 6, \"ContentDate\": {\"Start\": \"2024-12-21T17:34:36Z\", \"End\": \"2024-12-21T17:34:39Z\"}}, {\"Name\": \"s1a-wv1-ocn-vv-20241221t173450-20241221t173453-057093-0704a8-007.nc\", \"index\": 7, \"ContentDate\": {\"Start\": \"2024-12-21T17:34:50Z\", \"End\": \"2024-12-21T17:34:53Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t173505-20241221t173508-057093-0704a8-008.nc\", \"index\": 8, \"ContentDate\": {\"Start\": \"2024-12-21T17:35:05Z\", \"End\": \"2024-12-21T17:35:08Z\"}}]" + }, + { + "Name": "S1A_WV_OCN__2SSV_20241221T063740_20241221T063955_057086_070467_0F98.SAFE", + "Id": "232eff91-cdd5-4a9c-9d0c-c67a1ec4acd8", + "ContentDate": "{\"Start\": \"2024-12-21T06:37:40.503275Z\", \"End\": \"2024-12-21T06:39:55.264985Z\"}", + "GeoFootprint": "{\"type\": \"MultiPolygon\", \"coordinates\": [[[[178.404816, -29.766712], [178.618271, -29.71962], [178.568115, -29.541056], [178.355133, -29.588057], [178.404816, -29.766712]]], [[[179.994461, -28.467854], [180.0, -28.466545737174613], [180.0, -28.28197697389664], [179.943527, -28.295313], [179.994461, -28.467854]]], [[[-180.0, -28.466545737174613], [-179.801666, -28.419701], [-179.852936, -28.247248], [-180.0, -28.28197697389664], [-180.0, -28.466545737174613]]], [[[177.901077, -28.013386], [178.111328, -27.966906], [178.062653, -27.788027], [177.852859, -27.834421], [177.901077, -28.013386]]], [[[179.479889, -26.71768], [179.680542, -26.67037], [179.630737, -26.497107], [179.430405, -26.544333], [179.479889, -26.71768]]], [[[177.422653, -26.255707], [177.629288, -26.20993], [177.581985, -26.030941], [177.375778, -26.076635], [177.422653, -26.255707]]], [[[178.976395, -24.966898], [179.17424, -24.920357], [179.125977, -24.746952], [178.928421, -24.793413], [178.976395, -24.966898]]], [[[176.957108, -24.496502], [177.160431, -24.451385], [177.114319, -24.272079], [176.911392, -24.317118], [176.957108, -24.496502]]], [[[178.492538, -23.212112], [178.687683, -23.166328], [178.640778, -22.992655], [178.445892, -23.038364], [178.492538, -23.212112]]], [[[176.503189, -22.735729], [176.703476, -22.691233], [176.658493, -22.511806], [176.458557, -22.556231], [176.503189, -22.735729]]], [[[178.022034, -21.455568], [178.214706, -21.4105], [178.169067, -21.236683], [177.976639, -21.281679], [178.022034, -21.455568]]]]}", + "WVImagettes": "[{\"Name\": \"s1a-wv1-ocn-vv-20241221t063740-20241221t063743-057086-070467-001.nc\", \"index\": 1, \"ContentDate\": {\"Start\": \"2024-12-21T06:37:40Z\", \"End\": \"2024-12-21T06:37:43Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t063755-20241221t063758-057086-070467-002.nc\", \"index\": 2, \"ContentDate\": {\"Start\": \"2024-12-21T06:37:55Z\", \"End\": \"2024-12-21T06:37:58Z\"}}, {\"Name\": \"s1a-wv1-ocn-vv-20241221t063809-20241221t063812-057086-070467-003.nc\", \"index\": 3, \"ContentDate\": {\"Start\": \"2024-12-21T06:38:09Z\", \"End\": \"2024-12-21T06:38:12Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t063824-20241221t063827-057086-070467-004.nc\", \"index\": 4, \"ContentDate\": {\"Start\": \"2024-12-21T06:38:24Z\", \"End\": \"2024-12-21T06:38:27Z\"}}, {\"Name\": \"s1a-wv1-ocn-vv-20241221t063839-20241221t063842-057086-070467-005.nc\", \"index\": 5, \"ContentDate\": {\"Start\": \"2024-12-21T06:38:39Z\", \"End\": \"2024-12-21T06:38:42Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t063853-20241221t063856-057086-070467-006.nc\", \"index\": 6, \"ContentDate\": {\"Start\": \"2024-12-21T06:38:53Z\", \"End\": \"2024-12-21T06:38:56Z\"}}, {\"Name\": \"s1a-wv1-ocn-vv-20241221t063908-20241221t063911-057086-070467-007.nc\", \"index\": 7, \"ContentDate\": {\"Start\": \"2024-12-21T06:39:08Z\", \"End\": \"2024-12-21T06:39:11Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t063923-20241221t063925-057086-070467-008.nc\", \"index\": 8, \"ContentDate\": {\"Start\": \"2024-12-21T06:39:23Z\", \"End\": \"2024-12-21T06:39:25Z\"}}, {\"Name\": \"s1a-wv1-ocn-vv-20241221t063937-20241221t063940-057086-070467-009.nc\", \"index\": 9, \"ContentDate\": {\"Start\": \"2024-12-21T06:39:37Z\", \"End\": \"2024-12-21T06:39:40Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t063952-20241221t063955-057086-070467-010.nc\", \"index\": 10, \"ContentDate\": {\"Start\": \"2024-12-21T06:39:52Z\", \"End\": \"2024-12-21T06:39:55Z\"}}]" + }, + { + "Name": "S1A_WV_OCN__2SSV_20241221T024033_20241221T024051_057084_070452_58C0.SAFE", + "Id": "a8ab2e6b-9d16-4825-95c6-97cc956ca200", + "ContentDate": "{\"Start\": \"2024-12-21T02:40:33.868599Z\", \"End\": \"2024-12-21T02:40:51.453773Z\"}", + "GeoFootprint": "{\"type\": \"MultiPolygon\", \"coordinates\": [[[[47.039242, -4.94837], [47.077576, -4.767841], [46.893394, -4.727687], [46.854965, -4.908176], [47.039242, -4.94837]]], [[[45.204754, -5.469716], [45.242847, -5.294876], [45.062805, -5.255004], [45.024651, -5.429796], [45.204754, -5.469716]]]]}", + "WVImagettes": "[{\"Name\": \"s1a-wv1-ocn-vv-20241221t024033-20241221t024036-057084-070452-001.nc\", \"index\": 1, \"ContentDate\": {\"Start\": \"2024-12-21T02:40:33Z\", \"End\": \"2024-12-21T02:40:36Z\"}}, {\"Name\": \"s1a-wv2-ocn-vv-20241221t024048-20241221t024051-057084-070452-002.nc\", \"index\": 2, \"ContentDate\": {\"Start\": \"2024-12-21T02:40:48Z\", \"End\": \"2024-12-21T02:40:51Z\"}}]" + }, + { + "Name": "S1D_WV_OCN__2ASV_20260114T202903_20260114T205021_001035_000B05_82B5.SAFE", + "Id": "a2aa3f62-2402-41d9-8c85-94325b4812b1", + "ContentDate": "{\"Start\": \"2026-01-14T20:29:03.806293Z\", \"End\": \"2026-01-14T20:50:21.044384Z\"}", + "GeoFootprint": "{\"type\": \"MultiPolygon\", \"coordinates\": [[[[-48.662785, 64.343582], [-48.550373, 64.169281], [-48.142654, 64.21769], [-48.252518, 64.392227], [-48.662785, 64.343582]]], [[[-45.384644, 65.602875], [-45.29372, 65.432159], [-44.864391, 65.470512], [-44.952538, 65.641373], [-45.384644, 65.602875]]], [[[-49.857544, 66.055122], [-49.731819, 65.881416], [-49.296169, 65.932465], [-49.418995, 66.106438], [-49.857544, 66.055122]]], [[[-46.37751, 67.331528], [-46.275909, 67.161293], [-45.815773, 67.201401], [-45.914169, 67.371811], [-46.37751, 67.331528]]], [[[-51.223652, 67.755661], [-51.081676, 67.583153], [-50.612663, 67.637558], [-50.751324, 67.810387], [-51.223652, 67.755661]]], [[[-47.487995, 69.055489], [-47.373207, 68.886078], [-46.877827, 68.928261], [-46.988895, 69.0979], [-47.487995, 69.055489]]], [[[-52.705734, 69.453438], [-52.545334, 69.283768], [-52.043037, 69.341347], [-52.199696, 69.511398], [-52.705734, 69.453438]]], [[[-48.741261, 70.774765], [-48.609367, 70.605797], [-48.073444, 70.650436], [-48.20097, 70.819679], [-48.741261, 70.774765]]], [[[-54.43568, 71.137047], [-54.247955, 70.967026], [-53.705891, 71.028679], [-53.889271, 71.199158], [-54.43568, 71.137047]]], [[[-50.200787, 72.485855], [-50.046749, 72.317741], [-49.462769, 72.365463], [-49.611629, 72.533913], [-50.200787, 72.485855]]], [[[-56.451588, 72.804726], [-56.230755, 72.636536], [-55.643421, 72.703018], [-55.859257, 72.871758], [-56.451588, 72.804726]]], [[[-51.915756, 74.186943], [-51.732418, 74.020058], [-51.091572, 74.071571], [-51.268723, 74.238884], [-51.915756, 74.186943]]], [[[-58.846436, 74.452301], [-58.581791, 74.286423], [-57.942513, 74.358772], [-58.201412, 74.525322], [-58.846436, 74.452301]]], [[[-53.97551, 75.875549], [-53.752048, 75.710175], [-53.04311, 75.766426], [-53.259094, 75.932343], [-53.97551, 75.875549]]], [[[-61.736511, 76.074478], [-61.413395, 75.911995], [-60.715519, 75.991417], [-61.032127, 76.154732], [-61.736511, 76.074478]]], [[[-56.516022, 77.546852], [-56.235825, 77.383522], [-55.444733, 77.445793], [-55.715824, 77.609825], [-56.516022, 77.546852]]], [[[-65.264633, 77.663788], [-64.86235, 77.506279], [-64.101715, 77.593948], [-64.496864, 77.752487], [-65.264633, 77.663788]]], [[[-59.778828, 79.189911], [-59.415604, 79.029823], [-58.523754, 79.100105], [-58.875896, 79.261124], [-59.778828, 79.189911]]], [[[-69.889412, 79.183189], [-69.37043, 79.031685], [-68.534584, 79.131355], [-69.046249, 79.284172], [-69.889412, 79.183189]]], [[[-64.161163, 80.787613], [-63.672089, 80.633209], [-62.657024, 80.714577], [-63.133007, 80.870239], [-64.161163, 80.787613]]], [[[-75.795273, 80.631485], [-75.114304, 80.490181], [-74.213165, 80.603722], [-74.888275, 80.746681], [-75.795273, 80.631485]]], [[[-70.156715, 82.324028], [-69.464119, 82.177643], [-68.310364, 82.273689], [-68.988983, 82.421829], [-70.156715, 82.324028]]], [[[-83.709831, 81.950218], [-82.795464, 81.824272], [-81.856453, 81.955589], [-82.770447, 82.083542], [-83.709831, 81.950218]]], [[[-78.963097, 83.738403], [-77.932762, 83.606522], [-76.651566, 83.723396], [-77.673248, 83.857681], [-78.963097, 83.738403]]], [[[-94.403526, 83.060867], [-93.17672, 82.959442], [-92.285065, 83.112381], [-93.525406, 83.216026], [-94.403526, 83.060867]]], [[[-92.265045, 84.93367], [-90.706375, 84.830132], [-89.431458, 84.974724], [-91.008041, 85.08123], [-92.265045, 84.93367]]], [[[-108.385643, 83.867828], [-106.815392, 83.803772], [-106.158081, 83.977715], [-107.766136, 84.043602], [-108.385643, 83.867828]]], [[[-111.831566, 85.73053], [-109.622116, 85.677155], [-108.797379, 85.850441], [-111.088257, 85.906029], [-111.831566, 85.73053]]], [[[-125.241692, 84.23645], [-123.452209, 84.222359], [-123.272705, 84.408989], [-125.120407, 84.423546], [-125.241692, 84.23645]]], [[[-136.108994, 85.884323], [-133.691833, 85.902092], [-133.91124, 86.085129], [-136.437744, 86.066544], [-136.108994, 85.884323]]], [[[-142.717667, 84.081917], [-141.004883, 84.121307], [-141.390991, 84.304695], [-143.154434, 84.264069], [-142.717667, 84.081917]]], [[[-158.216827, 85.330719], [-156.311066, 85.411758], [-157.386871, 85.574829], [-159.339859, 85.490921], [-158.216827, 85.330719]]], [[[-158.049225, 83.441208], [-156.642105, 83.524971], [-157.417206, 83.691383], [-158.849457, 83.605492], [-158.049225, 83.441208]]], [[[-173.976425, 84.271698], [-172.696884, 84.391281], [-174.005249, 84.524582], [-175.285782, 84.402237], [-173.976425, 84.271698]]], [[[-170.040314, 82.436882], [-168.973236, 82.55143], [-169.899994, 82.696251], [-170.972565, 82.579529], [-170.040314, 82.436882]]], [[[175.607666, 82.93071], [176.45491, 83.070778], [175.21611, 83.179115], [174.381454, 83.036972], [175.607666, 82.93071]]], [[[-178.992813, 81.191544], [-178.198944, 81.325912], [-179.132584, 81.45079], [-179.922714, 81.314568], [-178.992813, 81.191544]]], [[[168.581039, 81.439133], [169.165207, 81.590401], [168.065155, 81.680443], [167.494888, 81.527679], [168.581039, 81.439133]]], [[[174.321671, 79.794518], [174.917709, 79.941216], [174.037659, 80.049561], [173.448486, 79.901382], [174.321671, 79.794518]]], [[[163.62236, 79.862747], [164.044128, 80.020164], [163.078201, 80.097023], [162.668594, 79.938515], [163.62236, 79.862747]]], [[[169.236465, 78.29895], [169.695023, 78.453796], [168.885788, 78.548988], [168.434601, 78.392967], [169.236465, 78.29895]]], [[[159.948257, 78.236198], [160.267487, 78.397881], [159.415512, 78.465034], [159.106384, 78.302536], [159.948257, 78.236198]]], [[[165.273636, 76.738335], [165.635925, 76.898849], [164.897217, 76.983635], [164.541824, 76.822182], [165.273636, 76.738335]]], [[[157.116898, 76.577156], [157.367081, 76.741409], [156.608932, 76.80127], [156.36702, 76.636406], [157.116898, 76.577156]]], [[[162.117172, 75.131844], [162.410263, 75.296272], [161.735275, 75.372871], [161.448334, 75.207687], [162.117172, 75.131844]]], [[[154.860443, 74.896263], [155.062729, 75.062286], [154.381577, 75.116562], [154.186111, 74.95005], [154.860443, 74.896263]]], [[[159.528442, 73.494576], [159.770721, 73.661858], [159.152313, 73.731789], [158.915421, 73.563889], [159.528442, 73.494576]]], [[[153.01059, 73.200134], [153.178391, 73.367393], [152.560852, 73.41729], [152.398727, 73.249649], [153.01059, 73.200134]]], [[[157.368042, 71.834084], [157.572021, 72.00338], [157.002899, 72.067871], [156.803589, 71.898071], [157.368042, 71.834084]]], [[[151.462357, 71.491653], [151.604584, 71.659752], [151.04007, 71.7062], [150.902588, 71.537804], [151.462357, 71.491653]]], [[[155.530716, 70.156097], [155.705597, 70.32711], [155.179382, 70.387108], [155.008591, 70.215675], [155.530716, 70.156097]]], [[[150.1306, 69.775864], [150.253754, 69.945099], [149.7341, 69.988739], [149.614975, 69.81926], [150.1306, 69.775864]]], [[[153.945389, 68.464485], [154.097305, 68.636604], [153.608414, 68.692871], [153.460052, 68.520401], [153.945389, 68.464485]]], [[[148.971786, 68.052773], [149.079956, 68.222672], [148.59845, 68.264053], [148.493759, 68.093948], [148.971786, 68.052773]]], [[[152.569275, 66.760864], [152.703262, 66.934074], [152.246262, 66.987328], [152.115417, 66.81382], [152.569275, 66.760864]]], [[[147.946472, 66.324829], [148.042694, 66.49514], [147.594055, 66.53466], [147.500824, 66.364189], [147.946472, 66.324829]]], [[[151.330826, 65.050705], [151.45018, 65.224754], [151.022415, 65.275269], [150.905792, 65.10096], [151.330826, 65.050705]]], [[[147.022522, 64.593407], [147.108673, 64.763145], [146.688965, 64.801094], [146.605408, 64.63121], [147.022522, 64.593407]]], [[[150.211441, 63.335754], [150.317108, 63.507744], [149.915894, 63.555794], [149.812607, 63.383591], [150.211441, 63.335754]]], [[[146.180954, 62.85733], [146.259735, 63.028347], [145.865265, 63.064976], [145.788757, 62.893848], [146.180954, 62.85733]]], [[[149.211441, 61.611404], [149.308624, 61.786091], [148.929977, 61.83218], [148.834946, 61.65731], [149.211441, 61.611404]]], [[[145.42775, 61.117565], [145.499893, 61.288181], [145.127106, 61.323822], [145.056961, 61.153111], [145.42775, 61.117565]]], [[[148.334213, 59.8778], [148.423706, 60.05352], [148.062836, 60.098324], [147.975281, 59.92244], [148.334213, 59.8778]]], [[[144.743271, 59.37389], [144.810257, 59.544937], [144.456116, 59.579845], [144.390915, 59.408722], [144.743271, 59.37389]]], [[[147.52179, 58.142101], [147.604614, 58.318443], [147.259918, 58.362103], [147.178848, 58.185619], [147.52179, 58.142101]]], [[[144.094894, 57.628525], [144.1577, 57.800514], [143.820541, 57.834766], [143.759354, 57.662708], [144.094894, 57.628525]]], [[[146.749207, 56.404427], [146.826157, 56.581097], [146.497345, 56.623554], [146.421982, 56.446762], [146.749207, 56.404427]]], [[[143.489288, 55.881184], [143.548325, 56.053448], [143.226456, 56.087162], [143.168869, 55.914845], [143.489288, 55.881184]]], [[[146.029099, 54.663364], [146.101059, 54.840431], [145.78656, 54.881847], [145.716034, 54.704678], [146.029099, 54.663364]]], [[[142.913635, 54.131813], [142.969254, 54.303791], [142.661255, 54.337048], [142.606918, 54.165028], [142.913635, 54.131813]]], [[[145.354294, 52.919586], [145.421951, 53.096973], [145.120407, 53.137497], [145.054047, 52.960014], [145.354294, 52.919586]]], [[[142.3685, 52.381042], [142.421509, 52.553734], [142.125992, 52.586628], [142.074158, 52.413902], [142.3685, 52.381042]]], [[[144.72583, 51.171177], [144.789749, 51.348751], [144.499481, 51.38858], [144.436752, 51.210926], [144.72583, 51.171177]]], [[[141.848068, 50.628967], [141.898666, 50.801815], [141.614502, 50.834423], [141.564987, 50.661549], [141.848068, 50.628967]]], [[[144.121582, 49.422588], [144.182281, 49.600544], [143.902832, 49.639683], [143.843216, 49.461658], [144.121582, 49.422588]]], [[[141.352097, 48.874962], [141.400497, 49.047668], [141.126572, 49.08007], [141.079132, 48.907345], [141.352097, 48.874962]]], [[[143.551147, 47.670593], [143.608856, 47.848328], [143.338959, 47.886925], [143.282242, 47.709126], [143.551147, 47.670593]]], [[[140.879654, 47.118885], [140.926331, 47.292019], [140.661591, 47.324303], [140.615799, 47.151157], [140.879654, 47.118885]]], [[[143.003006, 45.91713], [143.05835, 46.095535], [142.797409, 46.133629], [142.742981, 45.955166], [143.003006, 45.91713]]], [[[140.416351, 45.362778], [140.461426, 45.536076], [140.205338, 45.568241], [140.161072, 45.394936], [140.416351, 45.362778]]], [[[142.479324, 44.163372], [142.532349, 44.34164], [142.279617, 44.379299], [142.227417, 44.200985], [142.479324, 44.163372]]], [[[139.971802, 43.604782], [140.015533, 43.778419], [139.767197, 43.810543], [139.724197, 43.636902], [139.971802, 43.604782]]], [[[141.981491, 42.40443], [142.032547, 42.58276], [141.786774, 42.620155], [141.736481, 42.441784], [141.981491, 42.40443]]], [[[139.541809, 41.845501], [139.58429, 42.019184], [139.342972, 42.051319], [139.301163, 41.877636], [139.541809, 41.845501]]], [[[141.495041, 40.646774], [141.544479, 40.825703], [141.305542, 40.862801], [141.256805, 40.683838], [141.495041, 40.646774]]], [[[139.120316, 40.08551], [139.161743, 40.259411], [138.927032, 40.291565], [138.88623, 40.117668], [139.120316, 40.08551]]], [[[141.010635, 38.889427], [141.058395, 39.068214], [140.826691, 39.104885], [140.779572, 38.926067], [141.010635, 38.889427]]], [[[138.701904, 38.325748], [138.742065, 38.49855], [138.513687, 38.530704], [138.474075, 38.357906], [138.701904, 38.325748]]]]}", + "WVImagettes": null + } +] diff --git a/tests/test_sar_parquet.py b/tests/test_sar_parquet.py new file mode 100644 index 0000000..7c36007 --- /dev/null +++ b/tests/test_sar_parquet.py @@ -0,0 +1,228 @@ +"""CDSE daily parquet source: dateline merge, centres, logs of what is left out, a full run.""" + +import json +import logging +from pathlib import Path + +import numpy as np +import pyarrow as pa +import pyarrow.parquet as pq +import pytest + +from soba_matchups import sar_parquet as sp + +ROWS = json.loads((Path(__file__).parent / "data/parquet_rows.json").read_text()) +WEST_FIRST, EAST_FIRST, ORDINARY, A_CLASS = ROWS +ALL_MISSIONS = ("S1A", "S1B", "S1C", "S1D") + + +def km(a, b): + """Distance in km between two [lon, lat] points (small distances, longitude folded).""" + dx = ((a[0] - b[0] + 180) % 360 - 180) * np.cos(np.radians(a[1])) * 111.195 + return float(np.hypot(dx, (a[1] - b[1]) * 111.195)) + + +def write_day(folder, day, rows): + """Write a daily parquet file with the given rows (string columns, as the producer does).""" + folder.mkdir(exist_ok=True) + columns = {c: [r[c] for r in rows] for c in sp.COLUMNS} + pq.write_table(pa.table(columns), folder / f"s1_{day}.parquet") + + +def find(row): + """Return the day (YYYY-MM-DD) and the window around a row's SAFE.""" + day = json.loads(row["ContentDate"])["Start"][:10] + return day, f"{day}T00:00:00Z", f"{day}T23:59:59Z" + + +def test_west_half_first_and_east_half_first_are_both_merged(): + for row, pair in ((WEST_FIRST, (6, 7)), (EAST_FIRST, (1, 2))): + rings = sp._rings(row["GeoFootprint"]) + groups = sp.merge_dateline(rings) + assert len(groups) == len(json.loads(row["WVImagettes"])) == len(rings) - 1 + assert list(pair) in groups + assert all(len(g) == 1 for g in groups if list(pair) != g) + + +def test_halves_need_matching_cut_latitudes(): + rings = sp._rings(WEST_FIRST["GeoFootprint"]) + shifted = [[(x, y + 0.5) for x, y in rings[7]]] + assert not sp._are_halves(rings[6], shifted[0]) + assert sp._are_halves(rings[6], rings[7]) and sp._are_halves(rings[7], rings[6]) + assert not sp._are_halves(rings[0], rings[1]) + # A half that meets the cut at the right start but ends elsewhere is not the other half. + top = max(y for x, y in rings[7] if x == 180.0) + stretched = [(x, y + 0.5 if (x == 180.0 and y == top) else y) for x, y in rings[7]] + assert sp._cut_latitudes(stretched, 180)[0] == sp._cut_latitudes(rings[7], 180)[0] + assert not sp._are_halves(rings[6], stretched) + bottom = min(y for x, y in rings[7] if x == 180.0) + lowered = [(x, y - 0.5 if (x == 180.0 and y == bottom) else y) for x, y in rings[7]] + assert sp._cut_latitudes(lowered, 180)[-1] == sp._cut_latitudes(rings[7], 180)[-1] + assert not sp._are_halves(rings[6], lowered) + + +def test_centre_of_an_ordinary_polygon_is_the_mean_of_its_corners(): + ring = sp._rings(ORDINARY["GeoFootprint"])[0] + got = sp.centre([ring]) + assert km(got, [np.mean([x for x, _ in ring]), np.mean([y for _, y in ring])]) < 0.3 + + +def test_centre_of_a_dateline_imagette_lies_between_its_halves_in_range(): + rings = sp._rings(WEST_FIRST["GeoFootprint"]) + got = sp.centre([rings[6], rings[7]]) + assert -180 <= got[0] < 180 and abs(got[0]) > 179.8 + west, east = np.mean(rings[6], axis=0), np.mean(rings[7], axis=0) + assert abs(got[1] - (west[1] + east[1]) / 2) < 0.1 + # The same imagette in the other order gives the same centre. + assert km(got, sp.centre([rings[7], rings[6]])) < 0.01 + + +def test_six_vertex_polygon_uses_the_area_not_the_first_corners(): + # A square with two extra vertices on one side: the centroid stays at the square's centre. + ring = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.0, 1.0), (0.5, 1.0), (0.0, 1.0)] + assert km(sp.centre([ring]), [0.5, 0.5]) < 0.1 + + +def test_discovery_logs_safes_without_list_skipped_safes_and_unavailable_imagettes( + tmp_path, caplog +): + broken = dict(ORDINARY, Name=ORDINARY["Name"].replace("58C0", "AAAA")) + imagettes = json.loads(broken["WVImagettes"]) + broken["WVImagettes"] = json.dumps(imagettes[:1]) # two polygons, one imagette + day = find(ORDINARY)[0] + write_day(tmp_path / "parquet", day, [WEST_FIRST, ORDINARY, broken]) + archive = tmp_path / "esa" + archive.mkdir() + # Put one imagette file of ORDINARY in the archive; the others stay absent. + folder = sp._safe_folder(archive, ORDINARY["Name"]) / "measurement" + folder.mkdir(parents=True) + first = json.loads(ORDINARY["WVImagettes"])[0]["Name"] + (folder / first).write_bytes(b"") + with caplog.at_level(logging.WARNING, logger="soba_matchups.sar_parquet"): + products, prov = sp.discover_parquet( + f"{day}T00:00:00Z", f"{day}T23:59:59Z", ALL_MISSIONS, tmp_path / "parquet", archive + ) + assert {p["Name"] for p in products} == {WEST_FIRST["Name"], ORDINARY["Name"]} + assert [s["safe"] for s in prov["safes_skipped"]] == [broken["Name"]] + assert "imagettes" in prov["safes_skipped"][0]["reason"] + unavailable = {(u["safe"], u["index"]) for u in prov["imagettes_unavailable"]} + n_west = len(json.loads(WEST_FIRST["WVImagettes"])) + assert len(unavailable) == n_west + 1 # all of the dateline SAFE and the second ORDINARY file + assert (ORDINARY["Name"], 2) in unavailable and (ORDINARY["Name"], 1) not in unavailable + assert {u["reason"] for u in prov["imagettes_unavailable"]} == { + "SAFE not in the archive", + "file not in the archive", + } + ordinary = next(p for p in products if p["Name"] == ORDINARY["Name"]) + assert [i["available"] for i in ordinary["Imagettes"]] == [True, False] + text = caplog.text + assert "SAFEs skipped" in text and "imagettes whose file is not in the archive" in text + + +def test_a_safe_with_no_imagette_list_is_logged_and_left_out(tmp_path, caplog): + day = json.loads(A_CLASS["ContentDate"])["Start"][:10] + write_day(tmp_path / "parquet", day, [A_CLASS, ORDINARY]) + archive = tmp_path / "esa" + archive.mkdir() + with caplog.at_level(logging.WARNING, logger="soba_matchups.sar_parquet"): + products, prov = sp.discover_parquet( + f"{day}T00:00:00Z", f"{day}T23:59:59Z", ALL_MISSIONS, tmp_path / "parquet", archive + ) + assert products == [] # ORDINARY is from another day, A_CLASS has no list + assert prov["safes_without_imagettes"] == [ + {"safe": A_CLASS["Name"], "reason": "no imagette list in the parquet file"} + ] + assert "1 SAFEs without an imagette list" in caplog.text + + +def test_missing_columns_and_missing_files_are_reported(tmp_path): + day = find(ORDINARY)[0] + folder = tmp_path / "parquet" + folder.mkdir() + pq.write_table(pa.table({"Name": ["x"]}), folder / f"s1_{day}.parquet") + with pytest.raises(ValueError, match="lacks the columns"): + sp.discover_parquet(f"{day}T00:00:00Z", f"{day}T23:59:59Z", ALL_MISSIONS, folder, tmp_path) + with pytest.raises(ValueError, match="No daily parquet files"): + sp.discover_parquet( + "2000-01-02T00:00:00Z", "2000-01-02T01:00:00Z", ALL_MISSIONS, folder, tmp_path + ) + + +def test_a_changed_parquet_file_is_refused(tmp_path): + day = find(ORDINARY)[0] + write_day(tmp_path / "parquet", day, [ORDINARY]) + _, prov = sp.discover_parquet( + f"{day}T00:00:00Z", f"{day}T23:59:59Z", ALL_MISSIONS, tmp_path / "parquet", tmp_path + ) + sp.verify_parquet_inputs(prov) + write_day(tmp_path / "parquet", day, [WEST_FIRST]) + with pytest.raises(ValueError, match="changed"): + sp.verify_parquet_inputs(prov) + + +def point_ids(item): + """Return the sorted record ids of every point of an Item.""" + groups = item["properties"]["sarmatchup:matchups"] + return sorted(p["record_id"] for g in groups for p in g["points"]) + + +def synthetic_parquet(local_recipe): + """A daily file for the synthetic SAFE of conftest (its second imagette crosses 180).""" + from tests.conftest import SAFE + + root = local_recipe.parent + east = [[179.8, 10.0], [180.0, 10.0], [180.0, 11.0], [179.8, 11.0], [179.8, 10.0]] + west = [[-180.0, 10.0], [-179.8, 10.0], [-179.8, 11.0], [-180.0, 11.0], [-180.0, 10.0]] + square = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]] + names = sorted((root / "esa").rglob("*.nc")) + listed = [] + for index, path in enumerate(names, 1): + token = path.name.split("-") + start = f"2026-06-29T{token[4][9:11]}:{token[4][11:13]}:{token[4][13:15]}Z" + end = f"2026-06-29T{token[5][9:11]}:{token[5][11:13]}:{token[5][13:15]}Z" + listed.append( + {"Name": path.name, "index": index, "ContentDate": {"Start": start, "End": end}} + ) + row = { + "Name": SAFE, + "Id": "synthetic", + "ContentDate": json.dumps({"Start": "2026-06-29T12:00:00Z", "End": "2026-06-29T12:10:03Z"}), + "GeoFootprint": json.dumps( + {"type": "MultiPolygon", "coordinates": [[square], [east], [west]]} + ), + "WVImagettes": json.dumps(listed), + } + write_day(root / "parquet", "2026-06-29", [row]) + return root / "parquet" + + +def test_full_run_from_the_parquet_gives_the_same_matches_as_the_local_archive(local_recipe): + from soba_matchups.cli import main + + parquet = synthetic_parquet(local_recipe) + recipe = local_recipe.parent / "from_parquet.toml" + recipe.write_text( + local_recipe.read_text().replace('output_dir = "result"', 'output_dir = "from_parquet"') + + f'sar_parquet_root = "{parquet}"\n' + ) + assert main(["run", "--recipe", str(recipe)]) == 0 + root = local_recipe.parent + assert main(["run", "--recipe", str(local_recipe)]) == 0 + from_parquet = json.loads((root / "from_parquet/manifest.json").read_text()) + from_archive = json.loads((root / "result/manifest.json").read_text()) + assert from_parquet["discovery"]["source"] == "cdse_daily_parquet" + assert from_parquet["stac_items"] == from_archive["stac_items"] == 2 + assert ( + from_parquet["matchup_imagette_cell_rows"] + == from_archive["matchup_imagette_cell_rows"] + == 8 + ) + assert from_parquet["safes_without_imagettes"] == from_parquet["imagettes_unavailable"] == 0 + assert "daily parquet" in from_parquet["limitations"][0] + names = sorted(p.name for p in (root / "from_parquet/items").glob("*.json")) + assert names == sorted(p.name for p in (root / "result/items").glob("*.json")) + for name in names: + a = json.loads((root / "from_parquet/items" / name).read_text()) + b = json.loads((root / "result/items" / name).read_text()) + assert a["id"] == b["id"] # the same Item ids in both sources + assert point_ids(a) == point_ids(b) From 302ef2e4e145e9b74283ff6fc0dc8ba4f70499ef Mon Sep 17 00:00:00 2001 From: ilias <121384694+iliasreg@users.noreply.github.com> Date: Thu, 5 Nov 2026 14:14:43 +0100 Subject: [PATCH 2/3] Document sar_parquet_root in the docs --- docs/index.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/index.rst b/docs/index.rst index 393a758..d7e28df 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -33,6 +33,13 @@ with corresponding mission paths for S1B/C/D. Each extracted SAFE needs its ``stop_date`` is exclusive. ``output_dir`` must be new. Replace the previous ``sar_cache_dir`` / ``offline`` keys with ``sar_archive_root``. +To read the SAR metadata from the CDSE daily parquet files instead of opening the SAFEs, +add ``sar_parquet_root`` to the recipe: the folder of the ``s1_YYYY-MM-DD.parquet`` files, which +need the ``WVImagettes`` column. The footprint polygons give the imagette footprints and +centres. ``sar_archive_root`` is then used only to check which imagette files exist. The +manifest lists the SAFEs left out (no imagette list, or polygons and imagettes that disagree) +and the imagettes whose file is not in the archive, under ``discovery``. + Reference products ------------------ From 3b6f426dc8b786a4f8e83b2e1de3f43d83a473dd Mon Sep 17 00:00:00 2001 From: ilias <121384694+iliasreg@users.noreply.github.com> Date: Thu, 5 Nov 2026 16:19:00 +0100 Subject: [PATCH 3/3] Accept the archive, the parquet folder, or both as the SAR source A recipe needs sar_archive_root, sar_parquet_root, or both. The parquet folder alone skips the archive check and counts every imagette as available. --- docs/index.rst | 11 +++++--- src/soba_matchups/config.py | 20 ++++++++------ src/soba_matchups/sar_parquet.py | 19 +++++++++----- tests/test_config.py | 30 +++++++++++++++++++++ tests/test_sar_parquet.py | 45 ++++++++++++++++++++++++++++++++ 5 files changed, 107 insertions(+), 18 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index d7e28df..2761aaf 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -34,11 +34,14 @@ with corresponding mission paths for S1B/C/D. Each extracted SAFE needs its ``sar_cache_dir`` / ``offline`` keys with ``sar_archive_root``. To read the SAR metadata from the CDSE daily parquet files instead of opening the SAFEs, -add ``sar_parquet_root`` to the recipe: the folder of the ``s1_YYYY-MM-DD.parquet`` files, which +set ``sar_parquet_root`` in the recipe: the folder of the ``s1_YYYY-MM-DD.parquet`` files, which need the ``WVImagettes`` column. The footprint polygons give the imagette footprints and -centres. ``sar_archive_root`` is then used only to check which imagette files exist. The -manifest lists the SAFEs left out (no imagette list, or polygons and imagettes that disagree) -and the imagettes whose file is not in the archive, under ``discovery``. +centres. A recipe needs ``sar_archive_root``, ``sar_parquet_root``, or both. With the archive +alone the SAFEs are opened. With the parquet folder alone nothing is opened and every +imagette counts as available. With both, the metadata comes from the parquet files and the +archive is used only to check which imagette files exist. The manifest lists the SAFEs left +out (no imagette list, or polygons and imagettes that disagree) and the imagettes whose file is +not in the archive, under ``discovery``. Reference products ------------------ diff --git a/src/soba_matchups/config.py b/src/soba_matchups/config.py index 01a56f3..ebd0569 100644 --- a/src/soba_matchups/config.py +++ b/src/soba_matchups/config.py @@ -38,11 +38,11 @@ class Recipe: stop_date: date sar_missions: tuple[str, ...] time_tolerance_minutes: float - sar_archive_root: Path output_dir: Path instrument: str platform: str resolution: str + sar_archive_root: Path | None = None sar_parquet_root: Path | None = None def resolved(self): @@ -58,12 +58,16 @@ def load_recipe(path): with path.open("rb") as stream: data = tomllib.load(stream) # Reject unknown keys and missing required keys. - required = set(Recipe.__dataclass_fields__) - {"path", "sar_parquet_root"} - unknown = set(data) - required - {"sar_parquet_root"} + optional = {"sar_archive_root", "sar_parquet_root"} + required = set(Recipe.__dataclass_fields__) - {"path"} - optional + unknown = set(data) - required - optional if unknown or required - set(data): raise ValueError( f"Unknown recipe keys: {sorted(unknown)}; missing keys: {sorted(required - set(data))}" ) + # The SAR metadata comes from the archive, the daily parquet files, or both. + if not optional & set(data): + raise ValueError("Recipe needs sar_archive_root, sar_parquet_root, or both") product = (data["instrument"], data["platform"], data["resolution"]) if any(not isinstance(v, str) for v in product) or product not in REFERENCE_PRODUCTS: raise ValueError("Unsupported reference instrument / platform / resolution combination") @@ -100,11 +104,11 @@ def load_recipe(path): ): raise ValueError("sar_missions must contain distinct supported Sentinel-1 missions") # Make every path absolute, relative to the recipe file. - paths = ["reference_root", "sar_archive_root", "output_dir"] - inputs = ["reference_root", "sar_archive_root"] - if "sar_parquet_root" in data: - paths.append("sar_parquet_root") - inputs.append("sar_parquet_root") + paths = ["reference_root", "output_dir"] + inputs = ["reference_root"] + for key in sorted(optional & set(data)): + paths.append(key) + inputs.append(key) for key in paths: if not isinstance(data[key], str) or not data[key].strip(): raise ValueError(f"{key} must be a nonempty path") diff --git a/src/soba_matchups/sar_parquet.py b/src/soba_matchups/sar_parquet.py index 4b7753f..283394c 100644 --- a/src/soba_matchups/sar_parquet.py +++ b/src/soba_matchups/sar_parquet.py @@ -141,10 +141,10 @@ def _product(row, archive_root): if numbers != list(range(1, len(numbers) + 1)): raise ValueError("Imagette numbers are not 1..N") product_id = "local_" + hashlib.sha256(name.encode()).hexdigest() - folder = _safe_folder(archive_root, name) + folder = _safe_folder(archive_root, name) if archive_root is not None else None present = ( {p.name for p in (folder / "measurement").iterdir()} - if (folder / "measurement").is_dir() + if folder is not None and (folder / "measurement").is_dir() else None ) imagettes, unavailable = [], [] @@ -161,8 +161,11 @@ def _product(row, archive_root): "ContentDate": item["ContentDate"], "GeoFootprint": _geometry(member), "Centre": centre(member), - "href": (folder / relative).as_uri(), - "available": present is not None and item["Name"] in present, + # With no archive the file is not looked for: the link is the path inside the SAFE. + "href": (folder / relative).as_uri() + if folder is not None + else f"{name}/{relative}", + "available": folder is None or (present is not None and item["Name"] in present), } ) if not imagettes[-1]["available"]: @@ -181,7 +184,11 @@ def _product(row, archive_root): def discover_parquet( - start: str, end: str, missions: tuple[str, ...], parquet_root: Path, sar_archive_root: Path + start: str, + end: str, + missions: tuple[str, ...], + parquet_root: Path, + sar_archive_root: Path | None = None, ) -> tuple[list[dict], dict]: """Read the daily parquet files of the covered days and return the WV OCN products to match. @@ -200,7 +207,7 @@ def discover_parquet( "coverage": {"start": start, "end": end}, "missions": sorted(set(missions)), "sar_parquet_root": str(parquet_root), - "sar_archive_root": str(sar_archive_root), + "sar_archive_root": None if sar_archive_root is None else str(sar_archive_root), "coverage_scope": "products listed in the selected daily parquet files", "safes_without_imagettes": [], "safes_skipped": [], diff --git a/tests/test_config.py b/tests/test_config.py index a01f554..f66bdd7 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -85,3 +85,33 @@ def test_invalid_knmi_selection(local_recipe, old, new): local_recipe.write_text(local_recipe.read_text().replace(old, new)) with pytest.raises(ValueError, match="supported|reference_version"): load_recipe(local_recipe) + + +def test_a_recipe_needs_the_archive_the_parquet_folder_or_both(tmp_path): + from soba_matchups.config import load_recipe + + for name in ("reference", "esa", "parquet"): + (tmp_path / name).mkdir() + base = """instrument = "hscat" +platform = "hy-2b" +resolution = "25km" +reference_root = "reference" +start_date = "2020-01-10" +stop_date = "2020-01-11" +sar_missions = ["S1A"] +time_tolerance_minutes = 60 +output_dir = "run" +""" + recipe = tmp_path / "r.toml" + recipe.write_text(base) + with pytest.raises(ValueError, match="sar_archive_root, sar_parquet_root, or both"): + load_recipe(recipe) + recipe.write_text(base + 'sar_archive_root = "esa"\n') + config = load_recipe(recipe) + assert config.sar_archive_root == tmp_path / "esa" and config.sar_parquet_root is None + recipe.write_text(base + 'sar_parquet_root = "parquet"\n') + config = load_recipe(recipe) + assert config.sar_archive_root is None and config.sar_parquet_root == tmp_path / "parquet" + recipe.write_text(base + 'sar_archive_root = "esa"\nsar_parquet_root = "parquet"\n') + config = load_recipe(recipe) + assert config.sar_archive_root and config.sar_parquet_root diff --git a/tests/test_sar_parquet.py b/tests/test_sar_parquet.py index 7c36007..da10ada 100644 --- a/tests/test_sar_parquet.py +++ b/tests/test_sar_parquet.py @@ -226,3 +226,48 @@ def test_full_run_from_the_parquet_gives_the_same_matches_as_the_local_archive(l b = json.loads((root / "result/items" / name).read_text()) assert a["id"] == b["id"] # the same Item ids in both sources assert point_ids(a) == point_ids(b) + + +def test_without_an_archive_every_imagette_counts_as_available(tmp_path): + day = find(ORDINARY)[0] + write_day(tmp_path / "parquet", day, [WEST_FIRST, ORDINARY]) + products, prov = sp.discover_parquet( + f"{day}T00:00:00Z", f"{day}T23:59:59Z", ALL_MISSIONS, tmp_path / "parquet" + ) + assert {p["Name"] for p in products} == {WEST_FIRST["Name"], ORDINARY["Name"]} + assert prov["imagettes_unavailable"] == [] and prov["sar_archive_root"] is None + assert all(i["available"] for p in products for i in p["Imagettes"]) + + +def test_a_run_with_only_the_parquet_folder_matches_the_run_with_both(local_recipe): + from soba_matchups.cli import main + + parquet = synthetic_parquet(local_recipe) + root = local_recipe.parent + both = root / "both.toml" + both.write_text( + local_recipe.read_text().replace('output_dir = "result"', 'output_dir = "both"') + + f'sar_parquet_root = "{parquet}"\n' + ) + only = root / "only.toml" + only.write_text( + "\n".join( + line + for line in both.read_text() + .replace('output_dir = "both"', 'output_dir = "only"') + .splitlines() + if not line.startswith("sar_archive_root") + ) + + "\n" + ) + assert main(["run", "--recipe", str(both)]) == 0 + assert main(["run", "--recipe", str(only)]) == 0 + a = json.loads((root / "both/manifest.json").read_text()) + b = json.loads((root / "only/manifest.json").read_text()) + assert a["stac_items"] == b["stac_items"] == 2 + assert a["matchup_imagette_cell_rows"] == b["matchup_imagette_cell_rows"] == 8 + assert b["imagettes_unavailable"] == 0 and b["discovery"]["sar_archive_root"] is None + for item in sorted((root / "both/items").glob("*.json")): + x = json.loads(item.read_text()) + y = json.loads((root / "only/items" / item.name).read_text()) + assert x["id"] == y["id"] and point_ids(x) == point_ids(y)