From c8771449b607b9359db88ee010b104c541a648d0 Mon Sep 17 00:00:00 2001 From: Alessandro Amici Date: Sat, 17 Oct 2026 20:27:59 +0200 Subject: [PATCH 1/3] Name EOPF datatree nodes like the CPM Co-Authored-By: Claude Opus 5.5 --- docs/DATATREE.md | 28 +++++++++++++++++++--------- tests/test_30_eopf_product.py | 33 ++++++++++++++++++++++++++++++--- xarray_sentinel/eopf_product.py | 31 ++++++++++++++++++++++++++++--- 3 files changed, 77 insertions(+), 15 deletions(-) diff --git a/docs/DATATREE.md b/docs/DATATREE.md index 4526cdb..831afcc 100644 --- a/docs/DATATREE.md +++ b/docs/DATATREE.md @@ -84,7 +84,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 | |---|---| @@ -181,14 +183,22 @@ 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. calibration, noise_range: snake_case names, range coordinate from `pixel` as diff --git a/tests/test_30_eopf_product.py b/tests/test_30_eopf_product.py index d982522..ea62feb 100644 --- a/tests/test_30_eopf_product.py +++ b/tests/test_30_eopf_product.py @@ -17,6 +17,7 @@ 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" @@ -53,17 +54,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 +99,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"} diff --git a/xarray_sentinel/eopf_product.py b/xarray_sentinel/eopf_product.py index e1ebac3..f51f297 100644 --- a/xarray_sentinel/eopf_product.py +++ b/xarray_sentinel/eopf_product.py @@ -1,6 +1,8 @@ +import datetime import io import os import pathlib +import re from collections.abc import Callable from typing import Any @@ -19,6 +21,29 @@ "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 def split_bursts(pol_dataset: xr.Dataset) -> dict[str | None, xr.Dataset]: @@ -74,7 +99,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]] = {} @@ -111,9 +136,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, From 217bd24d8e15c364236ee04426dded0ae9036273 Mon Sep 17 00:00:00 2001 From: Alessandro Amici Date: Sat, 17 Oct 2026 21:19:36 +0200 Subject: [PATCH 2/3] Add EOPF antenna_pattern and coordinate_conversion in eopf_sentinel1 Co-Authored-By: Claude Opus 5.5 --- docs/DATATREE.md | 12 ++- tests/test_15_eopf_sentinel1.py | 57 +++++++++++++++ tests/test_30_eopf_product.py | 3 +- xarray_sentinel/eopf_product.py | 12 ++- xarray_sentinel/eopf_sentinel1.py | 66 +++++++++++++++++ xarray_sentinel/sentinel1.py | 118 +++--------------------------- 6 files changed, 155 insertions(+), 113 deletions(-) create mode 100644 tests/test_15_eopf_sentinel1.py create mode 100644 xarray_sentinel/eopf_sentinel1.py diff --git a/docs/DATATREE.md b/docs/DATATREE.md index 831afcc..0e10960 100644 --- a/docs/DATATREE.md +++ b/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 @@ -205,7 +210,12 @@ adds converters to the step 1 registry, steps 10 and 11 come last. `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. diff --git a/tests/test_15_eopf_sentinel1.py b/tests/test_15_eopf_sentinel1.py new file mode 100644 index 0000000..9d9f9e5 --- /dev/null +++ b/tests/test_15_eopf_sentinel1.py @@ -0,0 +1,57 @@ +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 diff --git a/tests/test_30_eopf_product.py b/tests/test_30_eopf_product.py index ea62feb..66f11f1 100644 --- a/tests/test_30_eopf_product.py +++ b/tests/test_30_eopf_product.py @@ -110,7 +110,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: @@ -124,7 +123,7 @@ 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), diff --git a/xarray_sentinel/eopf_product.py b/xarray_sentinel/eopf_product.py index f51f297..b0fdbb4 100644 --- a/xarray_sentinel/eopf_product.py +++ b/xarray_sentinel/eopf_product.py @@ -9,7 +9,7 @@ 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), @@ -46,6 +46,12 @@ def make_product_name(safe_name: str) -> str: 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, +} + + def split_bursts(pol_dataset: xr.Dataset) -> dict[str | None, xr.Dataset]: if "number_of_bursts" not in pol_dataset.attrs: return {None: pol_dataset} @@ -123,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( diff --git a/xarray_sentinel/eopf_sentinel1.py b/xarray_sentinel/eopf_sentinel1.py new file mode 100644 index 0000000..5f13f5a --- /dev/null +++ b/xarray_sentinel/eopf_sentinel1.py @@ -0,0 +1,66 @@ +from typing import Any + +import numpy as np +import xarray as xr + +from . import esa_safe + + +def open_antenna_pattern(annotation: dict[str, Any]) -> xr.Dataset: + patterns = esa_safe.get_tag( + annotation, "antennaPattern/antennaPatternList/antennaPattern" + ) + pattern_datasets = [] + for pattern in patterns: + slant_range_time = np.fromstring(pattern["slantRangeTime"]["$"], sep=" ") + elevation_angle = np.fromstring( + pattern["elevationAngle"]["$"], dtype=np.float32, sep=" " + ) + incidence_angle = np.fromstring( + pattern["incidenceAngle"]["$"], dtype=np.float32, sep=" " + ) + data_vars = { + "slant_range_time_ap": ("count", slant_range_time), + "elevation_angle": ("count", elevation_angle), + "incidence_angle": ("count", incidence_angle), + "swath": pattern["swath"], + "roll": pattern["roll"], + "terrain_height": pattern["terrainHeight"], + } + coords = {"count": np.arange(slant_range_time.size)} + pattern_dataset = xr.Dataset(data_vars=data_vars, coords=coords) + azimuth_time = [np.datetime64(pattern["azimuthTime"], "ns")] + pattern_datasets.append(pattern_dataset.expand_dims(azimuth_time=azimuth_time)) + return xr.concat(pattern_datasets, dim="azimuth_time", join="outer") + + +def open_coordinate_conversion(annotation: dict[str, Any]) -> xr.Dataset | None: + coordinate_conversion_list = esa_safe.get_tag( + annotation, "coordinateConversion/coordinateConversionList" + ) + if "coordinateConversion" not in coordinate_conversion_list: + return None + conversions = coordinate_conversion_list["coordinateConversion"] + srgr_coefficients = [ + np.fromstring(conversion["srgrCoefficients"]["$"], sep=" ") + for conversion in conversions + ] + grsr_coefficients = [ + np.fromstring(conversion["grsrCoefficients"]["$"], sep=" ") + for conversion in conversions + ] + data_vars = { + "gr0": ("azimuth_time", [c["gr0"] for c in conversions]), + "sr0": ("azimuth_time", [c["sr0"] for c in conversions]), + "slant_range_time": ( + "azimuth_time", + [c["slantRangeTime"] for c in conversions], + ), + "srgr_coefficients": (("azimuth_time", "degree"), srgr_coefficients), + "grsr_coefficients": (("azimuth_time", "degree"), grsr_coefficients), + } + coords = { + "azimuth_time": [np.datetime64(c["azimuthTime"], "ns") for c in conversions], + "degree": np.arange(srgr_coefficients[0].size), + } + return xr.Dataset(data_vars=data_vars, coords=coords) diff --git a/xarray_sentinel/sentinel1.py b/xarray_sentinel/sentinel1.py index 3ae226a..f00fdb8 100644 --- a/xarray_sentinel/sentinel1.py +++ b/xarray_sentinel/sentinel1.py @@ -24,7 +24,7 @@ import rasterio import xarray as xr -from . import conventions, esa_safe +from . import conventions, eopf_sentinel1, esa_safe SPEED_OF_LIGHT = 299_792_458 # m / s ONE_SECOND = np.timedelta64(1, "s") @@ -158,75 +158,6 @@ def open_reference_replica_dataset( return da -def open_antenna_pattern( - annotation_path: esa_safe.PathType, attrs: dict[str, Any] = {} -) -> xr.Dataset: - antenna_pattern_list = esa_safe.parse_tag_as_list( - annotation_path, "//antennaPattern/antennaPatternList/antennaPattern" - ) - - slant_range_time_list = [] - azimuth_time_list = [] - elevation_angle_list = [] - elevation_pattern_list = [] - incidence_angle_list = [] - terrain_height_list = [] - roll_list = [] - - for vector in antenna_pattern_list: - azimuth_time_list.append(vector["azimuthTime"]) - slant_range_time = np.fromstring( - vector["slantRangeTime"]["$"], dtype=np.float32, sep=" " - ) - slant_range_time_list.append(slant_range_time) - - elevation_angle = np.fromstring( - vector["elevationAngle"]["$"], dtype=np.float32, sep=" " - ) - elevation_angle_list.append(elevation_angle) - - elevation_pattern = np.fromstring( - vector["elevationPattern"]["$"], dtype=np.float32, sep=" " - ) - elevation_pattern_list.append( - elevation_pattern[::2] * 1j + elevation_pattern[1::2] - ) - - incidence_angle = np.fromstring( - vector["incidenceAngle"]["$"], dtype=np.float32, sep=" " - ) - incidence_angle_list.append(incidence_angle) - - terrain_height_list.append(vector["terrainHeight"]) - roll_list.append(vector["roll"]) - - slant_range_time_array = np.array(slant_range_time_list) - if (slant_range_time_array - slant_range_time_array[0]).any(): - raise ValueError( - "Unable to organise noise vectors in a regular line-pixel grid" - ) - data_vars = { - "elevationAngle": ( - ("azimuth_time", "slant_range_time"), - np.array(elevation_angle_list), - ), - # "elevationPattern": ( - # ( "azimuth_time", "slant_range_time"), np.array(elevation_pattern_list)), - "incidenceAngle": ( - ("azimuth_time", "slant_range_time"), - np.array(incidence_angle_list), - ), - "terrainHeight": ("azimuth_time", terrain_height_list), - "roll": ("azimuth_time", roll_list), - } - coords = { - "slant_range_time": slant_range_time_array[0], - "azimuth_time": [np.datetime64(dt, "ns") for dt in azimuth_time_list], - } - da = xr.Dataset(data_vars=data_vars, coords=coords, attrs=attrs) - return da - - def open_replica_dataset( annotation_path: esa_safe.PathType, attrs: dict[str, Any] = {} ) -> xr.Dataset: @@ -369,42 +300,15 @@ def open_noise_azimuth_dataset( def open_coordinate_conversion_dataset( annotation_path: esa_safe.PathType, attrs: dict[str, Any] = {} ) -> xr.Dataset: - coordinate_conversion = esa_safe.parse_tag_as_list( - annotation_path, "//coordinateConversionList/coordinateConversion" - ) - if len(coordinate_conversion) == 0: + annotation = esa_safe.parse_xml(annotation_path) + ds = eopf_sentinel1.open_coordinate_conversion(annotation) + if ds is None: raise TypeError("coordinateConversion tag not present in annotations") - - gr0 = [] - sr0 = [] - azimuth_time = [] - slant_range_time = [] - srgrCoefficients: list[list[float]] = [] - grsrCoefficients: list[list[float]] = [] - for values in coordinate_conversion: - sr0.append(values["sr0"]) - gr0.append(values["gr0"]) - azimuth_time.append(values["azimuthTime"]) - slant_range_time.append(values["slantRangeTime"]) - srgrCoefficients.append( - [float(v) for v in values["srgrCoefficients"]["$"].split()] - ) - grsrCoefficients.append( - [float(v) for v in values["grsrCoefficients"]["$"].split()] - ) - - coords: dict[str, Any] = {} - data_vars: dict[str, Any] = {} - coords["azimuth_time"] = [np.datetime64(dt, "ns") for dt in azimuth_time] - coords["degree"] = list(range(len(srgrCoefficients[0]))) - - data_vars["gr0"] = ("azimuth_time", gr0) - data_vars["sr0"] = ("azimuth_time", sr0) - data_vars["slant_range_time"] = ("azimuth_time", slant_range_time) - data_vars["srgrCoefficients"] = (("azimuth_time", "degree"), srgrCoefficients) - data_vars["grsrCoefficients"] = (("azimuth_time", "degree"), grsrCoefficients) - - return xr.Dataset(data_vars=data_vars, coords=coords, attrs=attrs) + ds = ds.rename( + srgr_coefficients="srgrCoefficients", grsr_coefficients="grsrCoefficients" + ) + ds.attrs.update(attrs) + return ds # taken from `np.cross` docs @@ -726,10 +630,7 @@ def find_available_groups( "gcp", # "replica", # "reference_replica", - # "antenna_pattern", ]: - if product_type == "GRD" and metadata_group == "antenna_pattern": - continue groups[f"{swath_pol_group}/{metadata_group}"] = [abspath] if product_type == "GRD": groups[f"{swath_pol_group}/coordinate_conversion"] = [abspath] @@ -1175,7 +1076,6 @@ def ground_range_to_slant_range_time( "noise_azimuth": open_noise_azimuth_dataset, "replica": open_replica_dataset, "reference_replica": open_reference_replica_dataset, - "antenna_pattern": open_antenna_pattern, } From 8dbc29cd7d3fc7819aeee028011d5d0cf8675625 Mon Sep 17 00:00:00 2001 From: Alessandro Amici Date: Sat, 17 Oct 2026 21:38:58 +0200 Subject: [PATCH 3/3] Keep line and pixel in the EOPF gcp dataset Co-Authored-By: Claude Opus 5.5 --- docs/DATATREE.md | 9 ++++-- tests/test_15_eopf_sentinel1.py | 9 ++++++ tests/test_30_eopf_product.py | 8 +++-- xarray_sentinel/eopf_product.py | 2 +- xarray_sentinel/eopf_sentinel1.py | 35 ++++++++++++++++++++++ xarray_sentinel/sentinel1.py | 49 +++++-------------------------- 6 files changed, 64 insertions(+), 48 deletions(-) diff --git a/docs/DATATREE.md b/docs/DATATREE.md index 0e10960..ae81bd9 100644 --- a/docs/DATATREE.md +++ b/docs/DATATREE.md @@ -167,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 @@ -205,7 +206,9 @@ adds converters to the step 1 registry, steps 10 and 11 come last. 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. diff --git a/tests/test_15_eopf_sentinel1.py b/tests/test_15_eopf_sentinel1.py index 9d9f9e5..e311b6a 100644 --- a/tests/test_15_eopf_sentinel1.py +++ b/tests/test_15_eopf_sentinel1.py @@ -55,3 +55,12 @@ def test_open_coordinate_conversion() -> None: ) 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() diff --git a/tests/test_30_eopf_product.py b/tests/test_30_eopf_product.py index 66f11f1..04763b6 100644 --- a/tests/test_30_eopf_product.py +++ b/tests/test_30_eopf_product.py @@ -24,7 +24,7 @@ / "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": [ @@ -36,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, {}) @@ -126,7 +130,7 @@ def test_open_datatree_groups( "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), diff --git a/xarray_sentinel/eopf_product.py b/xarray_sentinel/eopf_product.py index b0fdbb4..42c0aa4 100644 --- a/xarray_sentinel/eopf_product.py +++ b/xarray_sentinel/eopf_product.py @@ -14,7 +14,6 @@ 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), @@ -49,6 +48,7 @@ def make_product_name(safe_name: str) -> str: 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, } diff --git a/xarray_sentinel/eopf_sentinel1.py b/xarray_sentinel/eopf_sentinel1.py index 5f13f5a..28f737c 100644 --- a/xarray_sentinel/eopf_sentinel1.py +++ b/xarray_sentinel/eopf_sentinel1.py @@ -5,6 +5,16 @@ from . import esa_safe +GCP_VARIABLES = { + "latitude": "latitude", + "longitude": "longitude", + "height": "height", + "incidence_angle": "incidenceAngle", + "elevation_angle": "elevationAngle", + "azimuth_time_gcp": "azimuthTime", + "slant_range_time_gcp": "slantRangeTime", +} + def open_antenna_pattern(annotation: dict[str, Any]) -> xr.Dataset: patterns = esa_safe.get_tag( @@ -34,6 +44,31 @@ def open_antenna_pattern(annotation: dict[str, Any]) -> xr.Dataset: return xr.concat(pattern_datasets, dim="azimuth_time", join="outer") +def open_gcp(annotation: dict[str, Any]) -> xr.Dataset: + points = esa_safe.get_tag( + annotation, "geolocationGrid/geolocationGridPointList/geolocationGridPoint" + ) + lines = sorted({point["line"] for point in points}) + pixels = sorted({point["pixel"] for point in points}) + shape = (len(lines), len(pixels)) + grids = {name: np.full(shape, np.nan) for name in GCP_VARIABLES} + grids["azimuth_time_gcp"] = np.full(shape, np.datetime64("NaT", "ns")) + for point in points: + index = lines.index(point["line"]), pixels.index(point["pixel"]) + for name, tag in GCP_VARIABLES.items(): + grids[name][index] = point[tag] + # the annotated times of the first column and row, see quirk 1 in docs/DATATREE.md + coords = { + "azimuth_time": grids["azimuth_time_gcp"][:, 0], + "slant_range_time": grids["slant_range_time_gcp"][0], + "line": ("azimuth_time", lines), + "pixel": ("slant_range_time", pixels), + } + dims = ("azimuth_time", "slant_range_time") + data_vars = {name: (dims, grid) for name, grid in grids.items()} + return xr.Dataset(data_vars=data_vars, coords=coords) + + def open_coordinate_conversion(annotation: dict[str, Any]) -> xr.Dataset | None: coordinate_conversion_list = esa_safe.get_tag( annotation, "coordinateConversion/coordinateConversionList" diff --git a/xarray_sentinel/sentinel1.py b/xarray_sentinel/sentinel1.py index f00fdb8..37b2cf1 100644 --- a/xarray_sentinel/sentinel1.py +++ b/xarray_sentinel/sentinel1.py @@ -324,48 +324,13 @@ def is_clockwise(poly: list[tuple[float, float]]) -> bool: def open_gcp_dataset( annotation: esa_safe.PathOrFileType, attrs: dict[str, Any] = {} ) -> xr.Dataset: - geolocation_grid_points = esa_safe.parse_tag_as_list( - annotation, "//geolocationGridPoint" - ) - - azimuth_time = [] - slant_range_time = [] - line_set = set() - pixel_set = set() - for ggp in geolocation_grid_points: - if ggp["line"] not in line_set: - azimuth_time.append(np.datetime64(ggp["azimuthTime"], "ns")) - line_set.add(ggp["line"]) - if ggp["pixel"] not in pixel_set: - slant_range_time.append(ggp["slantRangeTime"]) - pixel_set.add(ggp["pixel"]) - shape = (len(azimuth_time), len(slant_range_time)) - dims = ("azimuth_time", "slant_range_time") - data_vars = { - "latitude": (dims, np.full(shape, np.nan), attrs), - "longitude": (dims, np.full(shape, np.nan), attrs), - "height": (dims, np.full(shape, np.nan), attrs), - "incidenceAngle": (dims, np.full(shape, np.nan), attrs), - "elevationAngle": (dims, np.full(shape, np.nan), attrs), - } - line = sorted(line_set) - pixel = sorted(pixel_set) - for ggp in geolocation_grid_points: - for var in data_vars: - j = line.index(ggp["line"]) - i = pixel.index(ggp["pixel"]) - data_vars[var][1][j, i] = ggp[var] - - ds = xr.Dataset( - data_vars=data_vars, - coords={ - "azimuth_time": azimuth_time, - "slant_range_time": slant_range_time, - "line": ("azimuth_time", line), - "pixel": ("slant_range_time", pixel), - }, - attrs=attrs, - ) + xml_dict = esa_safe.parse_xml(annotation) + ds = eopf_sentinel1.open_gcp(xml_dict) + ds = ds.drop_vars(["azimuth_time_gcp", "slant_range_time_gcp"]) + ds = ds.rename(incidence_angle="incidenceAngle", elevation_angle="elevationAngle") + for data_var in ds.data_vars.values(): + data_var.attrs.update(attrs) + ds.attrs.update(attrs) footprint = get_footprint_linestring(ds.azimuth_time, ds.slant_range_time, ds) geospatial_attrs = make_geospatial_attributes(footprint)