diff --git a/docs/DATATREE.md b/docs/DATATREE.md index 7677bec..3ee3a15 100644 --- a/docs/DATATREE.md +++ b/docs/DATATREE.md @@ -81,7 +81,9 @@ Other product types, from the CPM code and not yet checked on samples: - SM SLC: one node, the measurement is not split. - WV SLC: one node per vignette. - GRD: one node, range dimension `ground_range`, measurement `grd` uint16, plus - `conditions/coordinate_conversion`. + `conditions/coordinate_conversion`. We add `_FillValue` 0, the no-data value of the + border, so xarray decodes `grd` to float32 with the border as NaN, EOPF has the Zarr + `fill_value` 0 but no `_FillValue`. Chunks, from the CPM mapping: TOPS SLC `{azimuth_time: whole burst, slant_range_time: 4096}`, other types 4096 x 4096. diff --git a/tests/test_30_eopf_product.py b/tests/test_30_eopf_product.py index d99bf0d..1630aff 100644 --- a/tests/test_30_eopf_product.py +++ b/tests/test_30_eopf_product.py @@ -3,6 +3,7 @@ import warnings from typing import Any +import numpy as np import pytest import xarray as xr @@ -194,6 +195,25 @@ def test_open_datatree_grd() -> None: assert "ground_range" in node[group].dims +def test_open_datatree_grd_encoding(tmp_path: pathlib.Path) -> None: + res = eopf_product.open_datatree(GRD_IW_2, check_files_exist=True) + node = res["S01SIWGRD_20211223T051122_0025_B172_5371_039993"] + measurements = node["measurements"].to_dataset() + grd = measurements.grd.isel(azimuth_time=slice(2), ground_range=slice(3)).load() + assert grd.dtype == "uint16" + grd[0, 0, 0] = 65535 + + grd.to_dataset().to_zarr(tmp_path / "grd.zarr", consolidated=False) + + raw = xr.open_zarr(tmp_path / "grd.zarr", consolidated=False, mask_and_scale=False) + assert raw.grd.dtype == "uint16" + np.testing.assert_array_equal(raw.grd.values, grd.values) + decoded = xr.open_zarr(tmp_path / "grd.zarr", consolidated=False) + assert decoded.grd.dtype == "float32" + assert decoded.grd.values[0, 0, 0] == 65535 + xr.testing.assert_equal(decoded.grd, grd.astype("float32").where(grd != 0)) + + @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) diff --git a/xarray_sentinel/eopf_product.py b/xarray_sentinel/eopf_product.py index 4343d4f..5826821 100644 --- a/xarray_sentinel/eopf_product.py +++ b/xarray_sentinel/eopf_product.py @@ -123,6 +123,10 @@ def make_product_name(safe_name: str, common_attrs: dict[str, Any]) -> str: "conditions/gcp": eopf_sentinel1.open_gcp, "conditions/orbit": eopf_sentinel1.open_orbit, } +# the measurements are opened without mask_and_scale, so the GRD digital numbers stay +# uint16 from the TIFF to the Zarr as in EOPF. The fill value lets xarray decode them to +# float32 as rioxarray does, 0 is the no-data value of the GRD border, as the EOPF Zarr fill_value. +MEASUREMENT_ENCODINGS = {"grd": {"_FillValue": 0}} def split_measurements( @@ -141,6 +145,7 @@ def split_measurements( for burst_id, burst in bursts.items(): measurement = burst.measurement.rename(name) measurement.attrs = {"long_name": measurement.attrs["long_name"]} + measurement.encoding.update(MEASUREMENT_ENCODINGS.get(name, {})) # the CPM chunks, whole TOPS bursts along azimuth chunks = {"azimuth_time": azimuth_chunk, measurement.dims[1]: 4096} measurements[burst_id] = measurement.chunk(chunks).to_dataset() @@ -233,6 +238,7 @@ def open_datatree( xml_files[annotation_path], fs=fs, rasterio_chunks=rasterio_chunks, + mask_and_scale=False, ) closers.append(pol_dataset.close) pol_dataset = conventions.update_attributes(pol_dataset) diff --git a/xarray_sentinel/sentinel1.py b/xarray_sentinel/sentinel1.py index c5dbfbc..422bacf 100644 --- a/xarray_sentinel/sentinel1.py +++ b/xarray_sentinel/sentinel1.py @@ -498,14 +498,13 @@ def find_available_groups( def open_rasterio_dataarray( measurement: esa_safe.PathOrFileType, fs: fsspec.AbstractFileSystem | None, - chunks: dict[str, int] | None, + **kwargs: Any, ) -> xr.DataArray: + kwargs = {"lock": False} | kwargs # fsspec needs rasterio >= 1.3.0, but we allow earlier rasterio versions for local files if fs is None or isinstance(fs, fsspec.implementations.local.LocalFileSystem): try: - arr = xr.open_dataarray( - measurement, engine="rasterio", chunks=chunks, lock=False - ) + arr = xr.open_dataarray(measurement, engine="rasterio", **kwargs) except rasterio.RasterioIOError as ex: if "No such file" in str(ex): raise FileNotFoundError(str(ex)) from ex @@ -515,9 +514,7 @@ def open_rasterio_dataarray( # full fsspec support via `fs.open(measurement)` is now broken and we fall back # to hope that rasterio remote protocols works with fsspec urlpaths maybe_rasterio_urlpath = fs.unstrip_protocol(measurement) - arr = xr.open_dataarray( - maybe_rasterio_urlpath, engine="rasterio", chunks=chunks, lock=False - ) + arr = xr.open_dataarray(maybe_rasterio_urlpath, engine="rasterio", **kwargs) return arr @@ -542,6 +539,7 @@ def open_pol_dataset( attrs: dict[str, Any] = {}, gcp: xr.Dataset | None = None, rasterio_chunks: dict[str, int] | None = None, + mask_and_scale: bool = True, ) -> xr.Dataset: product_information = esa_safe.parse_tag(annotation, "//productInformation") image_information = esa_safe.parse_tag(annotation, "//imageInformation") @@ -648,7 +646,9 @@ def open_pol_dataset( else: raise ValueError(f"unknown projection {product_information['projection']}") - arr = open_rasterio_dataarray(measurement, fs, rasterio_chunks) + arr = open_rasterio_dataarray( + measurement, fs, chunks=rasterio_chunks, mask_and_scale=mask_and_scale + ) close = arr.close preferred_chunks = arr.encoding["preferred_chunks"]