default fill_value not masked when read from file

[kmuehlbauer]

What happened?

When reading a netcdf file wich has been created with fill_value=None (default) those data is not masked. If one is writing back to disk this manifests.

What did you expect to happen?

Values should be masked.

There seems to be a simple solution:

On read apply the netcdf default fill_value in the variables attributes before decoding if no _FillValue attribute is set. After decoding we could change that to np.nan for floating point types.

Minimal Complete Verifiable Example

import numpy as np
import netCDF4 as nc
import xarray as xr

with nc.Dataset("test-no-missing-01.nc", mode="w") as ds:
    x = ds.createDimension("x", 5)
    test = ds.createVariable("test", "f4", ("x",), fill_value=None)
    test[:4] = np.array([0.0, np.nan, 1.0, 8.0], dtype="f4")
with nc.Dataset("test-no-missing-01.nc") as ds:
    print(ds["test"])
    print(ds["test"][:])

with xr.open_dataset("test-no-missing-01.nc").load() as roundtrip:
    print(roundtrip)
    print(roundtrip["test"].attrs)
    print(roundtrip["test"].encoding)
    roundtrip.to_netcdf("test-no-missing-02.nc")
with nc.Dataset("test-no-missing-02.nc") as ds:
    print(ds["test"])
    print(ds["test"][:])

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

<class 'netCDF4._netCDF4.Variable'>
float32 test(x)
unlimited dimensions: 
current shape = (5,)
filling on, default _FillValue of 9.969209968386869e+36 used
[0.0 nan 1.0 8.0 --]

<xarray.Dataset>
Dimensions:  (x: 5)
Dimensions without coordinates: x
Data variables:
    test     (x) float32 0.0 nan 1.0 8.0 9.969e+36
{}
{'zlib': False, 'szip': False, 'zstd': False, 'bzip2': False, 'blosc': False, 'shuffle': False, 'complevel': 0, 'fletcher32': False, 'contiguous': True, 'chunksizes': None, 'source': 'test-no-missing-01.nc', 'original_shape': (5,), 'dtype': dtype('float32')}
<class 'netCDF4._netCDF4.Variable'>
float32 test(x)
    _FillValue: nan
unlimited dimensions: 
current shape = (5,)
filling on
[0.0 -- 1.0 8.0 9.969209968386869e+36]

Anything else we need to know?

The issue is similar to #7722 but is more intricate, as now the status of certain data values change from masked to some netcdf specific default value.

This is when only parts of the source dataset have been written to. Then the default fill_value get's delivered to the user but it is not backed by an _FillValue attribute.

Environment

INSTALLED VERSIONS

commit: None
python: 3.11.0 | packaged by conda-forge | (main, Jan 14 2023, 12:27:40) [GCC 11.3.0]
python-bits: 64
OS: Linux
OS-release: 5.14.21-150400.24.55-default
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: de_DE.UTF-8
LOCALE: ('de_DE', 'UTF-8')
libhdf5: 1.14.0
libnetcdf: 4.9.2

xarray: 2023.3.0
pandas: 1.5.3
numpy: 1.24.2
scipy: 1.10.1
netCDF4: 1.6.3
pydap: None
h5netcdf: 1.1.0
h5py: 3.8.0
Nio: None
zarr: 2.14.2
cftime: 1.6.2
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: None
dask: 2023.3.1
distributed: 2023.3.1
matplotlib: None
cartopy: None
seaborn: None
numbagg: None
fsspec: 2023.3.0
cupy: 11.6.0
pint: 0.20.1
sparse: None
flox: None
numpy_groupies: None
setuptools: 67.6.0
pip: 23.0.1
conda: None
pytest: 7.2.2
mypy: None
IPython: 8.11.0
sphinx: None

[dcherian]

See #5680 (comment)

To follow up, from a practical perspective, there are two problems with assuming that there are always "truly missing values" (case 2):

It makes it impossible to represent the full range of values in a data type, e.g., 255 for uint8 now means "missing".
Due to unfortunately limited options for representing missing data in NumPy, Xarray represents truly missing values in its data model with "NaN". This is more or less OK for floating point data, but means that integer data gets converted into floats. For example, uint8 would now get automatically converted into float32.

Both of these issues are problematic for faithful "round tripping" of Xarray data into netCDF and back. For this reason, Xarray needs an unambiguous way to know if a netCDF variable could contain semantically missing values. So far, we've used the presence of missing_value and _FillValue attributes for that.

[kmuehlbauer]

@dcherian Great, a duplicate. :-( Sorry I must have overlooked that one.

It's somewhat counter-intuitive to get differing results when using netcdf4-python and xarray. Would be a good idea to document this behaviour.

It looks like it might at least be resolved for floating point source data:

Let's take the above simple example. We have np.nan written to the file, but the netcdf representation on disk uses a default (undeclared by attribute) _FillValue for unwritten parts.

For the netcdf4-python user the np.nan will not be masked, but the unfilled parts will be masked.

For xarray the default fillvalue won't be masked, appearing as valid data, which it is not. On subsequent writes np.nan will be introduced as the new fillvalue (by attribute), effectively changing the meaning of the default fillvalues.

Wouldn't it make sense then, to transform these default fill values to np.nan on read too, instead of giving the a seemingly meaningful value? Maybe yet another keyword switch, use_default_fillvalues?

There should be at least a warning on read, in these situations, that there are undefined values in the dataset which were never written and which will not be masked.

If the dataset contains unwritten parts, and a default fillvalue is used, in turn meaning the data creator did this by purpose (by not setting a _FillValue) it can mean several things:

• The creators data does actually not have missing values which need declaring, but it means, that his data will get masked for default fillvalue entries (maybe they doesn't know about this, but that might be unlikely).
• The creator doesn't care at all, with same conclusion as above.
• The creator purposefully uses default fillvalue as missing value, since they use this as a means of saving disk space. But this could also be done, by just defining that as _FillValue attribute at creation time, if I`m not mistaken.

I'm still convinced this could be fixed for floating point data.

[dcherian]

Would be a good idea to document this behaviour.

+1

Maybe yet another keyword switch, use_default_fillvalues?

Adding mask_default_netcdf_fill_values: bool is probably a good idea.

I'm still convinced this could be fixed for floating point data.

Generally its worse if we obey some default fill values but not others, because it becomes quite confusing to a user.

[kmuehlbauer]

I'm still convinced this could be fixed for floating point data.

Generally its worse if we obey some default fill values but not others, because it becomes quite confusing to a user.

I think this depends from which side you look at it :-) My point here is, we do not have to submissively obey to default fill values, but just use them when decoding. This only need to happen if no _FillValue is attached to the variable. By doing this we ensure that these missing values are mapped to np.nan (as it is expected by users).

In further course we can just apply the xarray standard np.nan when writing out. We need to document that in that case exact roundtrip isn't possible (it also isn't currently possible, in this example).

Consider this example:

dtype = "f4"
with nc.Dataset("test-fillvalues-01.nc", mode="w") as ds:
    x = ds.createDimension("x", 10)
    test_fillval_fillon = ds.createVariable("test_fillval_fillon", dtype, ("x",), fill_value=nc.default_fillvals[dtype])
    test_fillval_fillon[:5] = np.array([0.0, nc.default_fillvals[dtype], np.nan, 1.0, 8.0], dtype=dtype)
    test_nofillval_fillon = ds.createVariable("test_nofillval_fillon", dtype, ("x",), fill_value=None)
    test_nofillval_fillon[:5] = np.array([0.0, nc.default_fillvals[dtype], np.nan, 1.0, 8.0], dtype=dtype)

with nc.Dataset("test-fillvalues-01.nc") as ds:
    print("\n read with netCDF4-python")
    print("---------------------------")
    print(ds["test_fillval_fillon"])
    print(ds["test_fillval_fillon"][:])
    print(ds["test_nofillval_fillon"])
    print(ds["test_nofillval_fillon"][:])
    
with xr.open_dataset("test-fillvalues-01.nc").load() as ds:
    print("\n read with xarray")
    print("---------------------------")
    print(ds["test_fillval_fillon"])
    print(ds["test_fillval_fillon"][:])
    print(ds["test_nofillval_fillon"])
    print(ds["test_nofillval_fillon"][:])

 read with netCDF4-python
---------------------------
<class 'netCDF4._netCDF4.Variable'>
float32 test_fillval_fillon(x)
    _FillValue: 9.96921e+36
unlimited dimensions: 
current shape = (10,)
filling on
[0.0 -- nan 1.0 8.0 -- -- -- -- --]
<class 'netCDF4._netCDF4.Variable'>
float32 test_nofillval_fillon(x)
unlimited dimensions: 
current shape = (10,)
filling on, default _FillValue of 9.969209968386869e+36 used
[0.0 -- nan 1.0 8.0 -- -- -- -- --]

 read with xarray-python
---------------------------
<xarray.DataArray 'test_fillval_fillon' (x: 10)>
array([ 0., nan, nan,  1.,  8., nan, nan, nan, nan, nan], dtype=float32)
Dimensions without coordinates: x
<xarray.DataArray 'test_nofillval_fillon' (x: 10)>
array([0.00000e+00, 9.96921e+36,         nan, 1.00000e+00, 8.00000e+00,
       9.96921e+36, 9.96921e+36, 9.96921e+36, 9.96921e+36, 9.96921e+36],
      dtype=float32)
Dimensions without coordinates: x

The only difference between these two variables is that on the first the _FillValue is declared, and on the other the default _FillValue is used. So if xarray obeys (by CF standard) the first it should also obey the second.

This might just work, if these cases the default fillvalue is used for decoding to np.nan, and declared that np.nan will be the new _FillValue. Does that make sense?

[kmuehlbauer]

Closing. Please follow up in #2742.