Repository: zarr-developers/VirtualiZarr

Issue:
TypeError: ufunc 'isnan' not supported for the input types
Trying to open a local copy of one of the files comprising the CWorthy OAE atlas (xref #132) (stored [in the cloud here](https://source.coop/repositories/cworthy/oae-efficiency-atlas/description)) with @sharkinsspatial 's HDF reader raises an error

```python
In [4]: vds = open_virtual_dataset('../experimentation/virtualizarr/oae/alk-forcing.000-1999-01.pop.h.347.nc', backend=HDFVirtualBackend)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 vds = open_virtual_dataset('../experimentation/virtualizarr/oae/alk-forcing.000-1999-01.pop.h.347.nc', backend=HDFVirtualBackend)

File ~/Documents/Work/Code/virtualizarr/virtualizarr/backend.py:217, in open_virtual_dataset(filepath, filetype, group, drop_variables, loadable_variables, decode_times, cftime_variables, indexes, virtual_array_class, virtual_backend_kwargs, reader_options, backend)
    214 if backend_cls is None:
    215     raise NotImplementedError(f"Unsupported file type: {filetype.name}")
--> 217 vds = backend_cls.open_virtual_dataset(
    218     filepath,
    219     group=group,
    220     drop_variables=drop_variables,
    221     loadable_variables=loadable_variables,
    222     decode_times=decode_times,
    223     indexes=indexes,
    224     virtual_backend_kwargs=virtual_backend_kwargs,
    225     reader_options=reader_options,
    226 )
    228 return vds

File ~/Documents/Work/Code/virtualizarr/virtualizarr/readers/hdf/hdf.py:64, in HDFVirtualBackend.open_virtual_dataset(filepath, group, drop_variables, loadable_variables, decode_times, indexes, virtual_backend_kwargs, reader_options)
     55 drop_variables, loadable_variables = check_for_collisions(
     56     drop_variables,
     57     loadable_variables,
     58 )
     60 filepath = validate_and_normalize_path_to_uri(
     61     filepath, fs_root=Path.cwd().as_uri()
     62 )
---> 64 virtual_vars = HDFVirtualBackend._virtual_vars_from_hdf(
     65     path=filepath,
     66     group=group,
     67     drop_variables=drop_variables + loadable_variables,
     68     reader_options=reader_options,
     69 )
     71 loadable_vars, indexes = open_loadable_vars_and_indexes(
     72     filepath,
     73     loadable_variables=loadable_variables,
   (...)
     78     decode_times=decode_times,
     79 )
     81 attrs = HDFVirtualBackend._get_group_attrs(
     82     path=filepath, reader_options=reader_options, group=group
     83 )

File ~/Documents/Work/Code/virtualizarr/virtualizarr/readers/hdf/hdf.py:354, in HDFVirtualBackend._virtual_vars_from_hdf(path, group, drop_variables, reader_options)
    352 if key not in drop_variables:
    353     if isinstance(g[key], Dataset):
--> 354         variable = HDFVirtualBackend._dataset_to_variable(path, g[key])
    355         if variable is not None:
    356             variables[key] = variable

File ~/Documents/Work/Code/virtualizarr/virtualizarr/readers/hdf/hdf.py:284, in HDFVirtualBackend._dataset_to_variable(path, dataset)
    282 if isinstance(fill_value, np.ndarray):
    283     fill_value = fill_value[0]
--> 284 if np.isnan(fill_value):
    285     fill_value = float("nan")
    286 if isinstance(fill_value, np.generic):

TypeError: ufunc 'isnan' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
```

The error message is not very helpful. Numpy could have at least told me what the type it received was, and the HDF reader could have added context about which variable it was trying to parse when it failed (so that I could choose to load it instead). 

Discussion:

Patch:
diff --git a/README.md b/README.md
index 0fc46b1..91f07fd 100644
--- a/README.md
+++ b/README.md
@@ -13,7 +13,6 @@
 [![PyPI - Downloads](https://img.shields.io/pypi/dm/virtualizarr?label=pypi%7Cdownloads)](https://pypistats.org/packages/virtualizarr)
 [![Conda - Downloads](https://img.shields.io/conda/d/conda-forge/virtualizarr
 )](https://anaconda.org/conda-forge/virtualizarr)
-[![slack](https://img.shields.io/badge/slack-virtualizarr-purple.svg?logo=slack)](https://earthmover-community.slack.com/archives/C08EXCE8ZQX)
 
 ## Cloud-Optimize your Scientific Data as Virtual Zarr stores, using xarray syntax.
 
diff --git a/ci/environment.yml b/ci/environment.yml
index c2a0be4..a114b65 100644
--- a/ci/environment.yml
+++ b/ci/environment.yml
@@ -15,7 +15,7 @@ dependencies:
   - universal_pathlib
   - hdf5plugin
   - numcodecs
-  - icechunk>=0.2.0
+  - icechunk>=0.1.1
   # Testing
   - codecov[toml]
   - pre-commit
@@ -37,4 +37,4 @@ dependencies:
   - astropy
   - pip
   - pip:
-      - imagecodecs-numcodecs==2024.6.1
+    - imagecodecs-numcodecs==2024.6.1
diff --git a/ci/upstream.yml b/ci/upstream.yml
index 31fbb2d..f312de2 100644
--- a/ci/upstream.yml
+++ b/ci/upstream.yml
@@ -28,7 +28,7 @@ dependencies:
   - pooch
   - fsspec
   - pip
-  - icechunk>=0.2.0
   - pip:
+      - git+https://github.com/earth-mover/icechunk.git@main#subdirectory=icechunk-python # Installs zarr-python v3.0.0 as dependency
       - git+https://github.com/fsspec/kerchunk.git@main
       - imagecodecs-numcodecs==2024.6.1
diff --git a/docs/contributing.md b/docs/contributing.md
index 3c41e26..45fc859 100644
--- a/docs/contributing.md
+++ b/docs/contributing.md
@@ -71,5 +71,5 @@ Anyone with commit privileges to the repository can issue a release, and you sho
 8. Edit the draft release notes for consistency.
 9. Select 'Publish' to publish the release. This should automatically upload the new release to [PyPI](https://pypi.org/project/virtualizarr/) and [conda-forge](https://anaconda.org/conda-forge/virtualizarr).
 10. Check that this has run successfully (PyPI should show the new version number very quickly, but conda-forge might take several hours).
-11. Create and merge a PR to add a new empty section to the `docs/releases.rst` for the next release in the future. See [this commit](https://github.com/zarr-developers/VirtualiZarr/commit/e3912f08e22f2e3230af6eb1a2aacb5728822fa1) for an example (you can assume the next release will be numbered `vX.Y.Z+1`, but the number doesn't actually matter).
+11. Create and merge a PR to add a new empty section to the `docs/releases.rst` for the next release in the future.
 12. (Optional) Advertise the release on social media 📣
diff --git a/docs/releases.rst b/docs/releases.rst
index 4c82080..4865764 100644
--- a/docs/releases.rst
+++ b/docs/releases.rst
@@ -1,38 +1,13 @@
 Release notes
 =============
-.. _v1.3.2:
-
-v1.3.2 (unreleased)
--------------------
-
-New Features
-~~~~~~~~~~~~
-
-Breaking changes
-~~~~~~~~~~~~~~~~
-
-Deprecations
-~~~~~~~~~~~~
-
-Bug fixes
-~~~~~~~~~
-
-Documentation
-~~~~~~~~~~~~~
-
-Internal Changes
-~~~~~~~~~~~~~~~~
-
-- Updates `store.set_virtual_ref` to `store.set_virtual_refs` in `write_manifest_virtual_refs` (:pull:`443`) By `Raphael Hagen <https://github.com/norlandrhagen>`_.
 
 .. _v1.3.1:
 
-v1.3.1 (18th Feb 2025)
-----------------------
+v1.3.1 (unreleased)
+-------------------
 
 New Features
 ~~~~~~~~~~~~
-- Examples use new Icechunk syntax
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
diff --git a/pyproject.toml b/pyproject.toml
index 8083856..e6256af 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -71,7 +71,7 @@ all_readers = [
 
 # writers
 icechunk = [
-    "icechunk>=0.2.0",
+    "icechunk>=0.1.1",
 ]
 kerchunk = ["fastparquet"]
 all_writers = [
diff --git a/virtualizarr/readers/hdf/hdf.py b/virtualizarr/readers/hdf/hdf.py
index 2b45cd9..fc55c1b 100644
--- a/virtualizarr/readers/hdf/hdf.py
+++ b/virtualizarr/readers/hdf/hdf.py
@@ -8,11 +8,13 @@ from typing import (
     List,
     Mapping,
     Optional,
+    Tuple,
     Union,
 )
 
 import numpy as np
 import xarray as xr
+from xarray.backends.zarr import FillValueCoder
 
 from virtualizarr.manifests import (
     ChunkEntry,
@@ -40,6 +42,20 @@ else:
     H5Dataset: Any = None
     H5Group: Any = None
 
+FillValueType = Union[
+    int,
+    float,
+    bool,
+    complex,
+    str,
+    np.integer,
+    np.floating,
+    np.bool_,
+    np.complexfloating,
+    bytes,  # For fixed-length string storage
+    Tuple[bytes, int],  # Structured type
+]
+
 
 class HDFVirtualBackend(VirtualBackend):
     @staticmethod
@@ -216,6 +232,29 @@ class HDFVirtualBackend(VirtualBackend):
 
         return [dim.removeprefix(group) for dim in dims]
 
+    @staticmethod
+    def _extract_cf_fill_value(
+        h5obj: Union[H5Dataset, H5Group],
+    ) -> Optional[FillValueType]:
+        """
+        Convert the _FillValue attribute from an HDF5 group or dataset into
+        encoding.
+
+        Parameters
+        ----------
+        h5obj : h5py.Group or h5py.Dataset
+            An h5py group or dataset.
+        """
+        fillvalue = None
+        for n, v in h5obj.attrs.items():
+            if n == "_FillValue":
+                if isinstance(v, np.ndarray) and v.size == 1:
+                    fillvalue = v.item()
+                else:
+                    fillvalue = v
+                fillvalue = FillValueCoder.encode(fillvalue, h5obj.dtype)  # type: ignore[arg-type]
+        return fillvalue
+
     @staticmethod
     def _extract_attrs(h5obj: Union[H5Dataset, H5Group]):
         """
@@ -240,14 +279,14 @@ class HDFVirtualBackend(VirtualBackend):
         for n, v in h5obj.attrs.items():
             if n in _HIDDEN_ATTRS:
                 continue
+            if n == "_FillValue":
+                continue
             # Fix some attribute values to avoid JSON encoding exceptions...
             if isinstance(v, bytes):
                 v = v.decode("utf-8") or " "
             elif isinstance(v, (np.ndarray, np.number, np.bool_)):
                 if v.dtype.kind == "S":
                     v = v.astype(str)
-                if n == "_FillValue":
-                    continue
                 elif v.size == 1:
                     v = v.flatten()[0]
                     if isinstance(v, (np.ndarray, np.number, np.bool_)):
@@ -258,7 +297,6 @@ class HDFVirtualBackend(VirtualBackend):
                 v = ""
             if v == "DIMENSION_SCALE":
                 continue
-
             attrs[n] = v
         return attrs
 
@@ -290,21 +328,19 @@ class HDFVirtualBackend(VirtualBackend):
         codecs = codecs_from_dataset(dataset)
         cfcodec = cfcodec_from_dataset(dataset)
         attrs = HDFVirtualBackend._extract_attrs(dataset)
+        cf_fill_value = HDFVirtualBackend._extract_cf_fill_value(dataset)
+        attrs.pop("_FillValue", None)
+
         if cfcodec:
             codecs.insert(0, cfcodec["codec"])
             dtype = cfcodec["target_dtype"]
             attrs.pop("scale_factor", None)
             attrs.pop("add_offset", None)
-            fill_value = cfcodec["codec"].decode(dataset.fillvalue)
         else:
             dtype = dataset.dtype
-            fill_value = dataset.fillvalue
-        if isinstance(fill_value, np.ndarray):
-            fill_value = fill_value[0]
-        if np.isnan(fill_value):
-            fill_value = float("nan")
-        if isinstance(fill_value, np.generic):
-            fill_value = fill_value.item()
+
+        fill_value = dataset.fillvalue.item()
+
         filters = [codec.get_config() for codec in codecs]
         zarray = ZArray(
             chunks=chunks,  # type: ignore
@@ -323,6 +359,8 @@ class HDFVirtualBackend(VirtualBackend):
             variable = xr.Variable(data=marray, dims=dims, attrs=attrs)
         else:
             variable = xr.Variable(data=np.empty(dataset.shape), dims=dims, attrs=attrs)
+        if cf_fill_value is not None:
+            variable.encoding["_FillValue"] = cf_fill_value
         return variable
 
     @staticmethod
diff --git a/virtualizarr/vendor/__init__.py b/virtualizarr/vendor/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/virtualizarr/writers/icechunk.py b/virtualizarr/writers/icechunk.py
index f61378a..b5ec09a 100644
--- a/virtualizarr/writers/icechunk.py
+++ b/virtualizarr/writers/icechunk.py
@@ -283,18 +283,17 @@ def generate_chunk_key(
     index: tuple[int, ...],
     append_axis: Optional[int] = None,
     existing_num_chunks: Optional[int] = None,
-) -> list[int]:
+) -> str:
     if append_axis and append_axis >= len(index):
         raise ValueError(
             f"append_axis {append_axis} is greater than the number of indices {len(index)}"
         )
-
-    return [
-        ind + existing_num_chunks
+    return "/".join(
+        str(ind + existing_num_chunks)
         if axis is append_axis and existing_num_chunks is not None
-        else ind
+        else str(ind)
         for axis, ind in enumerate(index)
-    ]
+    )
 
 
 def write_manifest_virtual_refs(
@@ -307,18 +306,13 @@ def write_manifest_virtual_refs(
     last_updated_at: Optional[datetime] = None,
 ) -> None:
     """Write all the virtual references for one array manifest at once."""
-    from icechunk import VirtualChunkSpec
 
-    if group.name == "/":
-        key_prefix = arr_name
-    else:
-        key_prefix = f"{group.name}/{arr_name}"
+    key_prefix = f"{group.name}/{arr_name}"
 
     # loop over every reference in the ChunkManifest for that array
     # TODO inefficient: this should be replaced with something that sets all (new) references for the array at once
     # but Icechunk need to expose a suitable API first
     # See https://github.com/earth-mover/icechunk/issues/401 for performance benchmark
-
     it = np.nditer(
         [manifest._paths, manifest._offsets, manifest._lengths],  # type: ignore[arg-type]
         flags=[
@@ -329,15 +323,17 @@ def write_manifest_virtual_refs(
         op_flags=[["readonly"]] * 3,  # type: ignore
     )
 
-    virtual_chunk_spec_list = [
-        VirtualChunkSpec(
-            index=generate_chunk_key(it.multi_index, append_axis, existing_num_chunks),
+    for path, offset, length in it:
+        # it.multi_index will be an iterator of the chunk shape
+        index = it.multi_index
+        chunk_key = generate_chunk_key(index, append_axis, existing_num_chunks)
+
+        # set each reference individually
+        store.set_virtual_ref(
+            # TODO it would be marginally neater if I could pass the group and name as separate args
+            key=f"{key_prefix}/c/{chunk_key}",  # should be of form 'group/arr_name/c/0/1/2', where c stands for chunks
             location=path.item(),
             offset=offset.item(),
             length=length.item(),
-            last_updated_at_checksum=last_updated_at,
+            checksum=last_updated_at,
         )
-        for path, offset, length in it
-    ]
-
-    store.set_virtual_refs(array_path=key_prefix, chunks=virtual_chunk_spec_list)


Test patch:
diff --git a/virtualizarr/tests/test_readers/conftest.py b/virtualizarr/tests/test_readers/conftest.py
index 8d0c199..4884db4 100644
--- a/virtualizarr/tests/test_readers/conftest.py
+++ b/virtualizarr/tests/test_readers/conftest.py
@@ -342,3 +342,64 @@ def non_coord_dim(tmpdir):
     ds = ds.drop_dims("dim3")
     ds.to_netcdf(filepath, engine="netcdf4")
     return filepath
+
+
+@pytest.fixture
+def scalar_fill_value_hdf5_file(tmpdir):
+    filepath = f"{tmpdir}/scalar_fill_value.nc"
+    f = h5py.File(filepath, "w")
+    data = np.random.randint(0, 10, size=(5))
+    fill_value = 42
+    f.create_dataset(name="data", data=data, chunks=True, fillvalue=fill_value)
+    return filepath
+
+
+compound_dtype = np.dtype(
+    [
+        ("id", "i4"),  # 4-byte integer
+        ("temperature", "f4"),  # 4-byte float
+    ]
+)
+
+compound_data = np.array(
+    [
+        (1, 98.6),
+        (2, 101.3),
+    ],
+    dtype=compound_dtype,
+)
+
+compound_fill = (-9999, -9999.0)
+
+fill_values = [
+    {"fill_value": -9999, "data": np.random.randint(0, 10, size=(5))},
+    {"fill_value": -9999.0, "data": np.random.random(5)},
+    {"fill_value": np.nan, "data": np.random.random(5)},
+    {"fill_value": False, "data": np.array([True, False, False, True, True])},
+    {"fill_value": "NaN", "data": np.array(["three"], dtype="S10")},
+    {"fill_value": compound_fill, "data": compound_data},
+]
+
+
+@pytest.fixture(params=fill_values)
+def cf_fill_value_hdf5_file(tmpdir, request):
+    filepath = f"{tmpdir}/cf_fill_value.nc"
+    f = h5py.File(filepath, "w")
+    dset = f.create_dataset(name="data", data=request.param["data"], chunks=True)
+    dim_scale = f.create_dataset(
+        name="dim_scale", data=request.param["data"], chunks=True
+    )
+    dim_scale.make_scale()
+    dset.dims[0].attach_scale(dim_scale)
+    dset.attrs["_FillValue"] = request.param["fill_value"]
+    return filepath
+
+
+@pytest.fixture
+def cf_array_fill_value_hdf5_file(tmpdir):
+    filepath = f"{tmpdir}/cf_array_fill_value.nc"
+    f = h5py.File(filepath, "w")
+    data = np.random.random(5)
+    dset = f.create_dataset(name="data", data=data, chunks=True)
+    dset.attrs["_FillValue"] = np.array([np.nan])
+    return filepath
diff --git a/virtualizarr/tests/test_readers/test_hdf/test_hdf.py b/virtualizarr/tests/test_readers/test_hdf/test_hdf.py
index 598848a..d5146e8 100644
--- a/virtualizarr/tests/test_readers/test_hdf/test_hdf.py
+++ b/virtualizarr/tests/test_readers/test_hdf/test_hdf.py
@@ -1,6 +1,7 @@
 from unittest.mock import patch
 
 import h5py  # type: ignore
+import numpy as np
 import pytest
 
 from virtualizarr import open_virtual_dataset
@@ -111,6 +112,36 @@ class TestDatasetToVariable:
         )
         assert var.attrs["attribute_name"] == "attribute_name"
 
+    def test_scalar_fill_value(self, scalar_fill_value_hdf5_file):
+        f = h5py.File(scalar_fill_value_hdf5_file)
+        ds = f["data"]
+        var = HDFVirtualBackend._dataset_to_variable(
+            scalar_fill_value_hdf5_file, ds, group=""
+        )
+        assert var.data.zarray.fill_value == 42
+
+    def test_cf_fill_value(self, cf_fill_value_hdf5_file):
+        f = h5py.File(cf_fill_value_hdf5_file)
+        ds = f["data"]
+        if ds.dtype.kind in "S":
+            pytest.xfail("Investigate fixed-length binary encoding in Zarr v3")
+        if ds.dtype.names:
+            pytest.xfail(
+                "To fix, structured dtype fill value encoding for Zarr backend"
+            )
+        var = HDFVirtualBackend._dataset_to_variable(
+            cf_fill_value_hdf5_file, ds, group=""
+        )
+        assert "_FillValue" in var.encoding
+
+    def test_cf_array_fill_value(self, cf_array_fill_value_hdf5_file):
+        f = h5py.File(cf_array_fill_value_hdf5_file)
+        ds = f["data"]
+        var = HDFVirtualBackend._dataset_to_variable(
+            cf_array_fill_value_hdf5_file, ds, group=""
+        )
+        assert not isinstance(var.encoding["_FillValue"], np.ndarray)
+
 
 @requires_hdf5plugin
 @requires_imagecodecs
diff --git a/virtualizarr/tests/test_readers/test_hdf/test_hdf_integration.py b/virtualizarr/tests/test_readers/test_hdf/test_hdf_integration.py
index 8583f21..bf1aceb 100644
--- a/virtualizarr/tests/test_readers/test_hdf/test_hdf_integration.py
+++ b/virtualizarr/tests/test_readers/test_hdf/test_hdf_integration.py
@@ -6,9 +6,11 @@ import virtualizarr
 from virtualizarr.readers.hdf import HDFVirtualBackend
 from virtualizarr.tests import (
     requires_hdf5plugin,
+    requires_icechunk,
     requires_imagecodecs,
     requires_kerchunk,
 )
+from virtualizarr.tests.test_integration import roundtrip_as_in_memory_icechunk
 
 
 @requires_kerchunk
@@ -53,8 +55,12 @@ class TestIntegration:
         vds.virtualize.to_kerchunk(kerchunk_file, format="json")
         roundtrip = xr.open_dataset(kerchunk_file, engine="kerchunk")
         xrt.assert_allclose(ds, roundtrip)
+        assert (
+            ds["temperature"].encoding["_FillValue"]
+            == roundtrip["temperature"].encoding["_FillValue"]
+        )
 
-    def test_non_coord_dim(self, tmpdir, non_coord_dim):
+    def test_non_coord_dim_roundtrip(self, tmpdir, non_coord_dim):
         ds = xr.open_dataset(non_coord_dim)
         vds = virtualizarr.open_virtual_dataset(
             non_coord_dim, backend=HDFVirtualBackend
@@ -63,3 +69,18 @@ class TestIntegration:
         vds.virtualize.to_kerchunk(kerchunk_file, format="json")
         roundtrip = xr.open_dataset(kerchunk_file, engine="kerchunk")
         xrt.assert_equal(ds, roundtrip)
+
+    @requires_icechunk
+    def test_cf_fill_value_roundtrip(self, tmpdir, cf_fill_value_hdf5_file):
+        ds = xr.open_dataset(cf_fill_value_hdf5_file, engine="h5netcdf")
+        if ds["data"].dtype in [float, object]:
+            pytest.xfail(
+                "To fix handle fixed-length and structured type fill value \
+                encoding in xarray zarr backend."
+            )
+        vds = virtualizarr.open_virtual_dataset(
+            cf_fill_value_hdf5_file,
+            backend=HDFVirtualBackend,
+        )
+        roundtrip = roundtrip_as_in_memory_icechunk(vds, tmpdir, decode_times=False)
+        xrt.assert_equal(ds, roundtrip)
diff --git a/virtualizarr/tests/test_writers/test_icechunk.py b/virtualizarr/tests/test_writers/test_icechunk.py
index 01a8750..b7600c2 100644
--- a/virtualizarr/tests/test_writers/test_icechunk.py
+++ b/virtualizarr/tests/test_writers/test_icechunk.py
@@ -399,9 +399,7 @@ def test_generate_chunk_key_no_offset():
     # Test case without any offset (append_axis and existing_num_chunks are None)
     index = (1, 2, 3)
     result = generate_chunk_key(index)
-    assert result == [1, 2, 3], (
-        "The chunk key should match the index without any offset."
-    )
+    assert result == "1/2/3", "The chunk key should match the index without any offset."
 
 
 def test_generate_chunk_key_with_offset():
@@ -412,7 +410,7 @@ def test_generate_chunk_key_with_offset():
     result = generate_chunk_key(
         index, append_axis=append_axis, existing_num_chunks=existing_num_chunks
     )
-    assert result == [1, 7, 3], "The chunk key should offset the second index by 5."
+    assert result == "1/7/3", "The chunk key should offset the second index by 5."
 
 
 def test_generate_chunk_key_zero_offset():
@@ -423,7 +421,7 @@ def test_generate_chunk_key_zero_offset():
     result = generate_chunk_key(
         index, append_axis=append_axis, existing_num_chunks=existing_num_chunks
     )
-    assert result == [4, 5, 6], (
+    assert result == "4/5/6", (
         "No offset should be applied when existing_num_chunks is 0."
     )