Build an EOPF-like datatree from a SAFE (gcp, antenna_pattern and coordinate_conversions)

docs/DATATREE.md

@@ -16,6 +16,11 @@ dt = xr.open_datatree("S1A_IW_SLC__1SDV_....SAFE", engine="sentinel-1")
 - Entry point: `Sentinel1Backend.open_datatree` and `open_groups_as_dict`.
 - Strict EOPF content: no groups that EOPF lacks, for now. See `azimuth_fm_rate` in
   [Open points](#open-points).
+- The EOPF datasets are the base: `eopf_sentinel1` builds them from the `parse_xml`
+  dicts with `get_tag`. The `sentinel1` openers become wrappers converting them back to
+  the xarray-sentinel layout, with one `parse_xml` per call until all are ported. The
+  xarray-sentinel layout may get the EOPF float32 where float64 doesn't matter. Openers
+  never exposed in the xarray-sentinel layout, like antenna pattern, are EOPF only.
 
 ## References
 
@@ -84,7 +89,9 @@ S0{sat}S{mode}{type}_{start}_{duration:04}_{platform}{N:03}_{unique_id}_{datatak
 ```
 
 All fields come from the SAFE name, except `burst_id` from `swathTiming` and
-`N = (absolute_orbit - 73) // 175 + 1`, see quirk 5. WV appends `_{WVn}_{image_number}`.
+`N = (absolute_orbit - 73) // 175 + 1`, see quirk 5. `mode` is the beam for SM, e.g.
+`S01SS3SLC`. `burst_id` is the burst index from 0 before IPF 3.40. WV appends
+`_{WVn}_{image_number}`.
 
 | SAFE | EOPF CPM 3.0 node |
 |---|---|
@@ -160,8 +167,9 @@ What we do in each case follows the decisions above.
 ## Open points
 
 - `line` and `pixel` coordinates: CPM 2.7 has them, the CPM 3.0 sample does not. The
-  CPM 3.0 code drops them from calibration and noise but builds them for gcp. Check on a
-  CPM 3.0 conversion, step 11, until then omit them.
+  CPM 3.0 code drops them from calibration and noise but builds them for gcp. We keep
+  them in gcp, as additional data is safe. For the other groups check on a CPM 3.0
+  conversion, step 11, until then omit them.
 - `azimuth_fm_rate`, needed by sarsen: in EOPF it is only
   `other_metadata[pol]["general_annotation"]["azimuth_fm_rate_list"]`.
 - `zarr_conventions`, `spatial:bbox` and `proj:code` on the CPM 3.0 sample come from the
@@ -181,21 +189,36 @@ What we do in each case follows the decisions above.
 One PR with tests per step. Step 1 comes first, steps 2 to 9 are independent and each
 adds converters to the step 1 registry, steps 10 and 11 come last.
 
-1. Skeleton and test harness. Rewrite `eopf_product.open_datatree` around a context per
-   swath and polarisation that parses the manifest and each XML once and opens each
-   measurement once. Build each group with a registry `{eopf_path: converter}`, one
-   node per swath and burst, stacking the per-polarisation results along
-   `polarization` except orbit and attitude. Move `parse_eopf_metadata` out of
-   `open_sentinel1_dataset`. Add a helper comparing node paths, variable names, dims
-   and dtypes, with fixtures of the CPM 3.0 sample structure.
-1. Node names: a function tested on the examples above.
+1. Done: skeleton and test harness.
+   - `eopf_product.open_datatree` reads each XML once and opens each measurement once.
+     It builds the groups from the `EOPF_GROUPS` registry
+     `{eopf_path: (group, opener)}`, one node per swath and burst, stacking
+     polarisations along `polarization` except orbit and attitude.
+   - `parse_eopf_metadata` is gone from `open_sentinel1_dataset`,
+     `make_sentinel1_stac_item` is now `eopf_metadata.build_stac_discovery`.
+   - `esa_safe.parse_xml` decodes a whole XML once, `esa_safe.get_tag` indexes it by
+     path. The annotation is decoded once for `build_other_metadata` and
+     `build_stac_discovery`, while the registry openers and `open_pol_dataset` still
+     parse their XML at every call. New converters take the `parse_xml` dict and use
+     `get_tag`.
+   - `tests/test_30_eopf_product.py` compares each group with the CPM 3.0 burst node in
+     `tests/data/eopf/`. Groups not converted yet are marked `PENDING`, each step
+     removes the marks of its groups.
+1. Done: node names, `eopf_product.make_product_name` tested on the examples above.
 1. orbit, attitude, doppler_centroid: from the existing openers, quirks 2 and 3.
-1. gcp: snake_case names, 2-D annotated times, quirk 1.
+1. Done: gcp, `eopf_sentinel1.open_gcp` with snake_case names, 2-D annotated times,
+   quirk 1, and `line` and `pixel` as additional coordinates. It is the base of
+   `sentinel1.open_gcp_dataset`.
 1. calibration, noise_range: snake_case names, range coordinate from `pixel` as
    `image_slant_range_time + pixel / range_sampling_rate` for SLC and
    `pixel * range_pixel_spacing` for GRD.
 1. noise_azimuth: from the existing opener, quirk 6.
-1. antenna_pattern, and GRD coordinate_conversion.
+1. Done: antenna_pattern, and GRD coordinate_conversion. `eopf_product.ANNOTATION_GROUPS`
+   holds the converters of the `parse_xml` annotation, a converter returns `None` when
+   the product lacks the group. The GRD antenna patterns of different swaths are padded
+   with NaN along `count`. `eopf_sentinel1.open_coordinate_conversion` is the base of
+   `sentinel1.open_coordinate_conversion_dataset`, `eopf_sentinel1.open_antenna_pattern`
+   is EOPF only.
 1. measurements: split TOPS SLC by burst with `crop_burst_dataset` on one lazy array per
    swath and polarisation, quirk 4. Rename to `slc` or `grd`, keep only `long_name`,
    set the CPM chunks. SM and GRD as one image, WV per vignette.

tests/test_15_eopf_sentinel1.py

@@ -0,0 +1,66 @@
+import pathlib
+
+import xarray as xr
+
+from xarray_sentinel import eopf_sentinel1, esa_safe
+
+DATA_FOLDER = pathlib.Path(__file__).parent / "data"
+
+SLC_IW1_VV_annotation = (
+    DATA_FOLDER
+    / "S1B_IW_SLC__1SDV_20210401T052622_20210401T052650_026269_032297_EFA4.SAFE"
+    / "annotation"
+    / "s1b-iw1-slc-vv-20210401t052624-20210401t052649-026269-032297-004.xml"
+)
+GRD_IW_VV_annotation = (
+    DATA_FOLDER
+    / "S1B_IW_GRDH_1SDV_20210401T052623_20210401T052648_026269_032297_ECC8.SAFE"
+    / "annotation"
+    / "s1b-iw-grd-vv-20210401t052623-20210401t052648-026269-032297-001.xml"
+)
+
+
+def test_open_antenna_pattern() -> None:
+    res = eopf_sentinel1.open_antenna_pattern(esa_safe.parse_xml(GRD_IW_VV_annotation))
+
+    assert dict(res.sizes) == {"azimuth_time": 27, "count": 679}
+    assert set(res.swath.values.tolist()) == {"IW1", "IW2", "IW3"}
+    assert res.slant_range_time_ap.count("count").values.min() == 563
+
+
+def test_open_coordinate_conversion() -> None:
+    # dims and dtypes from the CPM 3.0 S01SIWGRD mapping
+    expected = {
+        "azimuth_time": [["azimuth_time"], "datetime64[ns]"],
+        "degree": [["degree"], "int64"],
+        "gr0": [["azimuth_time"], "float64"],
+        "grsr_coefficients": [["azimuth_time", "degree"], "float64"],
+        "slant_range_time": [["azimuth_time"], "float64"],
+        "sr0": [["azimuth_time"], "float64"],
+        "srgr_coefficients": [["azimuth_time", "degree"], "float64"],
+    }
+
+    res = eopf_sentinel1.open_coordinate_conversion(
+        esa_safe.parse_xml(GRD_IW_VV_annotation)
+    )
+
+    assert isinstance(res, xr.Dataset)
+    assert {
+        str(name): [list(var.dims), str(var.dtype)]
+        for name, var in res.variables.items()
+    } == expected
+
+    res = eopf_sentinel1.open_coordinate_conversion(
+        esa_safe.parse_xml(SLC_IW1_VV_annotation)
+    )
+
+    assert res is None
+
+
+def test_open_gcp() -> None:
+    res = eopf_sentinel1.open_gcp(esa_safe.parse_xml(SLC_IW1_VV_annotation))
+
+    assert dict(res.sizes) == {"azimuth_time": 10, "slant_range_time": 21}
+    assert res.azimuth_time_gcp.dtype == "datetime64[ns]"
+    assert (res.azimuth_time == res.azimuth_time_gcp.isel(slant_range_time=0)).all()
+    assert (res.slant_range_time == res.slant_range_time_gcp.isel(azimuth_time=0)).all()

tests/test_30_eopf_product.py

@@ -17,13 +17,14 @@
     DATA_FOLDER
     / "S1A_S3_SLC__1SDV_20210401T152855_20210401T152914_037258_04638E_6001.SAFE"
 )
+SLC_IW_NODE = "S01SIWSLC_20210401T052622_0028_B150_EFA4_032297_IW1_0"
 EOPF_SLC_IW_NODE = (
     DATA_FOLDER
     / "eopf"
     / "S01SIWSLC_20260731T234807_0014_C050_C079_0116E8_IW3_258176.json"
 )
 
-# EOPF quirks we fix, see docs/DATATREE.md
+# EOPF quirks we fix and data we add, see docs/DATATREE.md
 EOPF_FIXES = {
     "conditions/doppler_centroid": {
         "fine_dce_azimuth_start_time": [
@@ -35,6 +36,10 @@
             "datetime64[ns]",
         ],
     },
+    "conditions/gcp": {
+        "line": [["azimuth_time"], "int64"],
+        "pixel": [["slant_range_time"], "int64"],
+    },
 }
 EXPECTED_SLC_IW_STRUCTURE = {
     group: variables | EOPF_FIXES.get(group, {})
@@ -53,17 +58,43 @@ def datatree_structure(dt: xr.DataTree) -> dict[str, dict[str, list[Any]]]:
     }
 
 
+@pytest.mark.parametrize(
+    "safe_name,expected",
+    [
+        [
+            "S1C_IW_SLC__1SDV_20260731T234807_20260731T234821_008794_0116E8_C079",
+            "S01SIWSLC_20260731T234807_0014_C050_C079_0116E8",
+        ],
+        [
+            "S1C_IW_SLC__1SDH_20260531T194042_20260531T194109_007902_01011D_EC26",
+            "S01SIWSLC_20260531T194042_0027_C045_EC26_01011D",
+        ],
+        [
+            "S1C_IW_GRDH_1SDV_20261017T111256_20261017T111323_009924_013C8F_ECF5",
+            "S01SIWGRD_20261017T111256_0027_C057_ECF5_013C8F",
+        ],
+    ],
+)
+def test_make_product_name(safe_name: str, expected: str) -> None:
+    assert eopf_product.make_product_name(safe_name) == expected
+
+
+def test_make_product_name_error() -> None:
+    with pytest.raises(ValueError):
+        eopf_product.make_product_name("S2A_MSIL1C_20210403T101021_N0300_R022_T33TUM")
+
+
 @pytest.fixture(scope="module")
 def slc_iw_node_structure() -> dict[str, dict[str, list[Any]]]:
     dt = eopf_product.open_datatree(SLC_IW, check_files_exist=True)
-    return datatree_structure(dt.children[f"{SLC_IW.stem}_IW1_0"])
+    return datatree_structure(dt.children[SLC_IW_NODE])
 
 
 def test_open_datatree(capsys: pytest.CaptureFixture[str]) -> None:
     res = eopf_product.open_datatree(SLC_S3, check_files_exist=True)
 
     assert isinstance(res, xr.DataTree)
-    assert list(res.children) == [f"{SLC_S3.stem}_S3"]
+    assert list(res.children) == ["S01SS3SLC_20210401T152855_0019_A213_6001_04638E"]
     assert capsys.readouterr().out == ""
 
 
@@ -72,7 +103,7 @@ def test_open_datatree_bursts() -> None:
 
     assert len(res.children) == 19
 
-    node = res[f"{SLC_IW.stem}_IW1_0"]
+    node = res[SLC_IW_NODE]
     slc = node["measurements/slc"]
 
     assert node.attrs["other_metadata"].keys() == {"VH", "VV"}
@@ -83,7 +114,6 @@ def test_open_datatree_bursts() -> None:
     assert "polarization" not in node["conditions/orbit"].dims
 
 
-@PENDING
 def test_open_datatree_groups(
     slc_iw_node_structure: dict[str, dict[str, list[Any]]],
 ) -> None:
@@ -97,10 +127,10 @@ def test_open_datatree_groups(
         "conditions",
         "quality",
         pytest.param("measurements", marks=PENDING),
-        pytest.param("conditions/antenna_pattern", marks=PENDING),
+        "conditions/antenna_pattern",
         pytest.param("conditions/attitude", marks=PENDING),
         pytest.param("conditions/doppler_centroid", marks=PENDING),
-        pytest.param("conditions/gcp", marks=PENDING),
+        "conditions/gcp",
         pytest.param("conditions/orbit", marks=PENDING),
         pytest.param("quality/calibration", marks=PENDING),
         pytest.param("quality/noise_azimuth", marks=PENDING),

xarray_sentinel/eopf_product.py

@@ -1,24 +1,55 @@
+import datetime
 import io
 import os
 import pathlib
+import re
 from collections.abc import Callable
 from typing import Any
 
 import fsspec
 import xarray as xr
 
-from . import conventions, eopf_metadata, esa_safe, sentinel1
+from . import conventions, eopf_metadata, eopf_sentinel1, esa_safe, sentinel1
 
 EOPF_GROUPS: dict[str, tuple[str, Callable[..., xr.Dataset]]] = {
     "conditions/attitude": ("attitude", sentinel1.open_attitude_dataset),
     "conditions/doppler_centroid": ("dc_estimate", sentinel1.open_dc_estimate_dataset),
-    "conditions/gcp": ("gcp", sentinel1.open_gcp_dataset),
     "conditions/orbit": ("orbit", sentinel1.open_orbit_dataset),
     "quality/calibration": ("calibration", sentinel1.open_calibration_dataset),
     "quality/noise_azimuth": ("noise_azimuth", sentinel1.open_noise_azimuth_dataset),
     "quality/noise_range": ("noise_range", sentinel1.open_noise_range_dataset),
 }
 POLARIZATION_INDEPENDENT_GROUPS = {"conditions/attitude", "conditions/orbit"}
+SAFE_NAME_RE = re.compile(
+    r"S(\d)(\w)_(\w{2})_(\w{3})\w_\w{4}_(\w{15})_(\w{15})_(\d{6})_(\w{6})_(\w{4})"
+)
+SAFE_TIME_FORMAT = "%Y%m%dT%H%M%S"
+
+
+def make_product_name(safe_name: str) -> str:
+    match = SAFE_NAME_RE.fullmatch(safe_name)
+    if match is None:
+        raise ValueError(f"cannot parse SAFE name {safe_name!r}")
+    sat, platform, mode, product_type, start, stop, orbit, datatake_id, unique_id = (
+        match.groups()
+    )
+    start_time = datetime.datetime.strptime(start, SAFE_TIME_FORMAT)
+    stop_time = datetime.datetime.strptime(stop, SAFE_TIME_FORMAT)
+    duration = (stop_time - start_time).seconds
+    # not the relative orbit, see quirk 5 in docs/DATATREE.md
+    cycle = (int(orbit) - 73) // 175 + 1
+    product_name = (
+        f"S0{sat}S{mode}{product_type}_{start}_{duration:04}_{platform}{cycle:03}"
+        f"_{unique_id}_{datatake_id}"
+    )
+    return product_name
+
+
+ANNOTATION_GROUPS: dict[str, Callable[[dict[str, Any]], xr.Dataset | None]] = {
+    "conditions/antenna_pattern": eopf_sentinel1.open_antenna_pattern,
+    "conditions/coordinate_conversion": eopf_sentinel1.open_coordinate_conversion,
+    "conditions/gcp": eopf_sentinel1.open_gcp,
+}
 
 
 def split_bursts(pol_dataset: xr.Dataset) -> dict[str | None, xr.Dataset]:
@@ -74,7 +105,7 @@ def open_datatree(
             swath, _, pol = group.partition("/")
             swath_pols.setdefault(swath, []).append(pol)
 
-    product_name = pathlib.PurePosixPath(product_path).stem
+    product_name = make_product_name(pathlib.PurePosixPath(product_path).stem)
     dt = xr.DataTree()
     for swath, pols in swath_pols.items():
         bursts: dict[str | None, dict[str, xr.Dataset]] = {}
@@ -98,6 +129,10 @@ def open_datatree(
                     pol_datasets.setdefault(eopf_group, {})[pol] = ds.drop_attrs(
                         deep=False
                     )
+            for eopf_group, converter in ANNOTATION_GROUPS.items():
+                converted = converter(annotations[pol])
+                if converted is not None:
+                    pol_datasets.setdefault(eopf_group, {})[pol] = converted
             other_metadata[pol] = eopf_metadata.build_other_metadata(annotations[pol])
 
         stac_discovery = eopf_metadata.build_stac_discovery(
@@ -111,9 +146,9 @@ def open_datatree(
                 node_groups[eopf_group] = stack_polarizations(datasets)
 
         for burst_id, burst_pols in bursts.items():
-            node = f"{product_name}_{swath}"
+            node = product_name
             if burst_id is not None:
-                node += f"_{burst_id}"
+                node += f"_{swath}_{burst_id}"
             dt[node] = xr.Dataset(
                 attrs={
                     "other_metadata": other_metadata,