FillValue issues with kerchunk 1e37=> 9.99999993e+36

[rsignell-usgs]

In my original NetCDF4 file (HDF5 under the hood) the fill value _FillValue is set to 1.e+37 for 32bit floats, and when loaded into xarray with decode_cf=True, the data values with 1e37 are correctly interpreted as fill values and set to NaN.

In the kerchunk JSON, the _FillValue in the JSON is set to 9.99999993e+36 and when loaded into xarray with decode_cf=True, the data values are not interpreted as fill values, and show up as 9.99999993e+36 (not 1e37 and not NaN).

Pretty confused on what is going on here!

Here is a Reproducible Notebook

cc @peterm790

[martindurant]

My hunch was, that this is due to writing the value as text in the attributes of zarr. But that raises the question of how this is normally handled outside of kerchunk, and why decode_cf matters.

[martindurant]

@rsignell-usgs , do you mind using xarray's to_zarr on that dataset to see what .zarray/.zattrs it creates? You can set the output chunk store (as opposed to the metadata store) to be a dict or even a custom dict that doesn't store anything (see below) to avoid writing things.

class NoOpDict(dict):
    def __setitem__(self, *_):
        pass

metadata_strore = {}
chunk_store = NoOpDict()

[keewis]

this may be #105: if you look at the NetCDF file with h5py, you can see that _FillValue is array([1.e+37], dtype=float32) while the fillvalue property has a value of 0.0

[rsignell-usgs]

@keewis, what is "fillvalue"? I've never heard of that attribute being used for fill values. That must be an HDF5-specific thing?

@martindurant , when I write Zarr, it's even more confusing. The fill_value is set to 9.999999e+36, but when read back into xarray with cf_decode=False, the value is 1e37. ????

See cells [9] and [11] here:
https://nbviewer.org/gist/2400089529be8bc3924aacfd3b2485e6

[keewis]

That must be an HDF5-specific thing?

yes, exactly. You can examine it with

In [1]: import fsspec
   ...: import h5py
   ...: 
   ...: fs = fsspec.filesystem(
   ...:     "s3", anon=True, client_kwargs={"endpoint_url": "https://mghp.osn.xsede.org"}
   ...: )
   ...: url = "s3://rsignellbucket1/COAWST/coawst_us_20220101_01.nc"
   ...: 
   ...: f = h5py.File(fs.open(url))
   ...: f["temp"].fillvalue
Out[1]: 0.0

For float values it is used to set zarr's fill_value attribute, which AFAIU the zarr engine will use to replace _FillValue.

You might have to synchronise fill_value with _FillValue in the references (e.g. by passing the hacky function from #105 (comment) to MultiZarrToZarr's preprocess). For example:

In [1]: import fsspec
   ...: 
   ...: fs = fsspec.filesystem(
   ...:     "s3", anon=True, client_kwargs={"endpoint_url": "https://mghp.osn.xsede.org"}
   ...: )
   ...: url = "s3://rsignellbucket1/COAWST/coawst_us_20220101_01.nc"
   ...: 
   ...: import kerchunk.hdf
   ...: 
   ...: hdf = kerchunk.hdf.SingleHdf5ToZarr(fs.open(url), url=url)
   ...: references = hdf.translate()
       
In [2]: import xarray as xr
   ...: 
   ...: r_opts = dict(anon=True, client_kwargs={"endpoint_url": "https://mghp.osn.xsede.org"})

In [3]: fs = fsspec.filesystem(
   ...:     "reference", fo=references, remote_protocol="s3", remote_options=r_opts
   ...: )
   ...: m = fs.get_mapper("")
   ...: 
   ...: ds = xr.open_dataset(
   ...:     m,
   ...:     engine="zarr",
   ...:     chunks={},
   ...:     backend_kwargs=dict(consolidated=False),
   ...:     decode_cf=False,
   ...: )

In [4]: ds.temp.attrs
Out[4]: 
{'_FillValue': 0.0,
 'coordinates': 'lon_rho lat_rho s_rho ocean_time',
 'field': 'temperature, scalar, series',
 'grid': 'grid',
 'location': 'face',
 'long_name': 'potential temperature',
 'time': 'ocean_time',
 'units': 'Celsius'}

In [5]: ds.temp[0, 0, 0, 0].values
Out[5]: array(1.e+37, dtype=float32)

In [6]: decoded = xr.decode_cf(ds)
   ...: decoded.temp[0, 0, 0, 0].values
Out[6]: array(1.e+37, dtype=float32)

In [7]: import ujson
   ...: import itertools
   ...: 
   ...: 
   ...: def correct_fill_values(data):
   ...:     def fix_variable(values):
   ...:         zattrs = values[".zattrs"]
   ...: 
   ...:         if "_FillValue" not in zattrs:
   ...:             return values
   ...: 
   ...:         _FillValue = zattrs["_FillValue"]
   ...:         if values[".zarray"]["fill_value"] != _FillValue:
   ...:             values[".zarray"]["fill_value"] = _FillValue
   ...: 
   ...:         return values
   ...: 
   ...:     refs = data["refs"]
   ...:     prepared = (
   ...:         (tuple(key.split("/")), value) for key, value in refs.items() if "/" in key
   ...:     )
   ...:     filtered = (
   ...:         (key, ujson.loads(value))
   ...:         for key, value in prepared
   ...:         if key[1] in (".zattrs", ".zarray")
   ...:     )
   ...:     key = lambda i: i[0][0]
   ...:     grouped = (
   ...:         (name, {n[1]: v for n, v in group})
   ...:         for name, group in itertools.groupby(sorted(filtered, key=key), key=key)
   ...:     )
   ...:     fixed = ((name, fix_variable(var)) for name, var in grouped)
   ...:     flattened = {
   ...:         f"{name}/{item}": ujson.dumps(data, indent=4)
   ...:         for name, var in fixed
   ...:         for item, data in var.items()
   ...:     }
   ...:     data["refs"] = dict(sorted((refs | flattened).items()))
   ...:     return data

In [8]: corrected_refs = correct_fill_values(references)

In [9]: fs = fsspec.filesystem(
   ...:     "reference", fo=corrected_refs, remote_protocol="s3", remote_options=r_opts
   ...: )
   ...: m = fs.get_mapper("")
   ...: 
   ...: ds = xr.open_dataset(
   ...:     m,
   ...:     engine="zarr",
   ...:     chunks={},
   ...:     backend_kwargs=dict(consolidated=False),
   ...:     decode_cf=False,
   ...: )
   ...: display(ds.temp.attrs, ds.temp[0, 0, 0, 0].values)
{'_FillValue': 1e+37,
 'coordinates': 'lon_rho lat_rho s_rho ocean_time',
 'field': 'temperature, scalar, series',
 'grid': 'grid',
 'location': 'face',
 'long_name': 'potential temperature',
 'time': 'ocean_time',
 'units': 'Celsius'}
array(1.e+37, dtype=float32)

In [10]: decoded = xr.decode_cf(ds)
    ...: decoded.temp[0, 0, 0, 0].values
Out[10]: array(nan, dtype=float32)

[keewis]

The fill_value is set to 9.999999e+36, but when read back into xarray with cf_decode=False, the value is 1e37

That might be a precision / formatting issue:

In [33]: np.float32(9.999999933815813e36)
Out[33]: 1e+37

In [34]: np.float64(9.999999933815813e36)
Out[34]: 9.999999933815813e+36

[peterm790]

I have set up a reproducible notebook here manipulating some fill_values: https://gist.github.com/peterm790/6e6ddb2df0f095056837d21c531bbe57

It seems that only fill_value in .zarray is considered by xarray when opening as zarr from a kerchunk reference file and _FillValue in .zattrs is not considered.

The fill_value in @rsignell-usgs example is 0 (which is incorrect) so I am not sure if manually setting it to 1e37 would correctly render Nan, but might be worth investigating. This example didn't reproduce the 1e37=> 9.99999993e+36 despite setting dtype = npfloat32

[keewis]

just for reference, the code at the end of #177 (comment) successfully changes the references to work as expected (i.e. you get nan values instead of the fill value with decode_cf=True)

[martindurant]

@keewis , you mean the correct_fill_values function?

[keewis]

yes, although the the implementation feels a bit clunky

[rsignell-usgs]

@rabernat, hate to drag you into everything, but I bet you will immediately know what to do here, as you've been involved with xarray, the CF conventions & the xarray zarr writer.

The whole issue is explained in this This Reproducible Notebook (without credentials needed, as it uses OSN!)

[rabernat]

Rather than using string formatting to assess equality, I would be looking at raw bytes, e.g.

actual_value = ds.temp[0,0,0,0].data
actual_value_bytes = actual_value.tobytes()
assert actual_value_bytes == b'\xc2\xbd\xf0|'

I dug into this and noticed that the kerchunk dataset is converting the fill value to a float64 type, while the original (working) zarr on disk is converting it to a float32 (same as original netCDF).

Digging into this reveals some pretty weird things about how netCDF, Zarr, and Xarray handle fill values. Note, for example, that the on-disk zarr does not actually keep the fill value in the attributes at all:

! cat foo.zarr/temp/.zattrs
{
    "_ARRAY_DIMENSIONS": [
        "ocean_time",
        "s_rho",
        "eta_rho",
        "xi_rho"
    ],
    "coordinates": "lon_rho lat_rho s_rho ocean_time",
    "field": "temperature, scalar, series",
    "grid": "grid",
    "location": "face",
    "long_name": "potential temperature",
    "time": "ocean_time",
    "units": "Celsius"
}

Compare this to both the Kerchunk Zarr and the original NetCDF. Instead, the fill value lives in the zarr array attributes.

! cat foo.zarr/temp/.zarray
{
    "chunks": [
        1,
        4,
        84,
        448
    ],
    "compressor": {
        "blocksize": 0,
        "clevel": 5,
        "cname": "lz4",
        "id": "blosc",
        "shuffle": 1
    },
    "dtype": "<f4",
    "fill_value": 9.999999933815813e+36,
    "filters": null,
    "order": "C",
    "shape": [
        2,
        16,
        336,
        896
    ],
    "zarr_format": 2
}

Thus a key problem here is that there are redundant / overlapping mechanisms for specifying fill values. When we implemented the Zarr backend, we decided to lift the Zarr array fill_value into the dataset attributes when opening a new zarr array:
https://github.com/pydata/xarray/blob/d7931f9014a26e712ff5f30c4082cf0261f045d3/xarray/backends/zarr.py#L453-L456
(You can also see how writing is handled here.)

I am no longer sure that this was the right choice. The zarr definition of fill_value, is defined as

fill_value:
A scalar value providing the default value to use for uninitialized portions of the array, or null if no fill_value is to be used.

It is only ever used by zarr if you request data from a chunk that has not been initialized. It is NOT a mask. This is not quite the same is the NUG FillValue:

The _FillValue attribute specifies the fill value used to pre-fill disk space allocated to the variable. Such pre-fill occurs unless nofill mode is set using nc_set_fill(). The fill value is returned when reading values that were never written. If _FillValue is defined then it should be scalar and of the same type as the variable. If the variable is packed using scale_factor and add_offset attributes (see below), the _FillValue attribute should have the data type of the packed data.

It is not necessary to define your own _FillValue attribute for a variable if the default fill value for the type of the variable is adequate. However, use of the default fill value for data type byte is not recommended. Note that if you change the value of this attribute, the changed value applies only to subsequent writes; previously written data are not changed.

Generic applications often need to write a value to represent undefined or missing values. The fill value provides an appropriate value for this purpose because it is normally outside the valid range and therefore treated as missing when read by generic applications. It is legal (but not recommended) for the fill value to be within the valid range.

The kerchunk Zarr now has default fill_value of 0.0 in .zarray

print(fs.cat('temp/.zarray').decode("utf-8"))
{
    "chunks": [
        1,
        8,
        168,
        448
    ],
    "compressor": {
        "id": "zlib",
        "level": 1
    },
    "dtype": "<f4",
    "fill_value": 0.0,
    "filters": [
        {
            "elementsize": 4,
            "id": "shuffle"
        }
    ],
    "order": "C",
    "shape": [
        12,
        16,
        336,
        896
    ],
    "zarr_format": 2
}

---

Okay so why does any of this matter? Because whether we use the zarr array internal fill_value or the attribute _FillValue seems to affect how these are decoded from json to binary data. Observe

z_array = zarr.open("foo.zarr")["temp"]
on_disk_zarr_fill_value = z_array.fill_value
assert z_array[0, 0, 0, 0] == on_disk_zarr_fill_value
on_disk_zarr_fill_value_bytes = on_disk_zarr_fill_value.tobytes()
print(type(on_disk_zarr_fill_value), on_disk_zarr_fill_value_bytes)
# -> <class 'numpy.float32'> b'\xc2\xbd\xf0|'

but

kerhunk_z_array = zarr.open(m)['temp']
kerchunk_attrs_fill_value = kerhunk_z_array.attrs["_FillValue"]
print(type(kerchunk_attrs_fill_value))  # -> float

However, I am mystified by why this still works:

assert kerhunk_z_array[0, 0, 0, 0] == kerchunk_attrs_fill_value

... but cf-decoding does not.

I am currently looking into the xarray coding code to understand what is going wrong.

In the meantime, one thing you might try is to modify the kerchunk dataset to use the zarr Array fill_value property, instead of the CF-style attribute.

[rabernat]

Ok something even more obvious.

dsk = xr.open_dataset(m, engine="zarr", chunks={}, 
                     backend_kwargs=dict(consolidated=False), decode_cf=False)
dsk.temp._FillValue  # -> 0.0

So the zarr array fill_value is actually overwriting the fill value already in the attributes by Xarray! Manually overriding does the trick

dsk.temp.attrs['_FillValue'] = kerchunk_attrs_fill_value
xr.decode_cf(dsk).temp[0, 0, 0, 0].values  # -> NaN

This may be an xarray bug. The problem is ultimately the duplicate and inconsistent fill_values.

[rabernat]

This is the specific place where the overwriting happens: https://github.com/pydata/xarray/blob/d7931f9014a26e712ff5f30c4082cf0261f045d3/xarray/backends/zarr.py#L455-L456

So there are two simultaneous paths forward:

• Consider this an xarray bug, and make a PR to xarray to do something when it discovers two inconsistent fill values attributes. The most likely action would be to raise an error or warning.
• Regardless of that, kerchunk should strive to not produce inconsistent data (different values of array .fill_value and attribute _FillValue.

[rsignell-usgs]

@rabernat you are awesome! Thanks for showing us the path(s)!