Use uint16 on disk and still open as a float32 unless mask_and_scale=…

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.

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)

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)

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"]