Numeric scalar variable attributes (including fill_value, scale_factor, add_offset) are 1-d instead of 0-d with h5netcdf engine, triggering ValueError: non-broadcastable output on application when loading single elements

[gerritholl]

What happened:

When I try to open a NetCDF file using the h5netcdf engine, accessing a single data point before scale factors have been applied results in ValueError: non-broadcastable output operand with shape () doesn't match the broadcast shape (1,). The MCVE (see below) results in:

Traceback (most recent call last):
  File "mwe93.py", line 4, in <module>
    ds["Rad"][400, 300].load()
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/core/dataarray.py", line 808, in load
    ds = self._to_temp_dataset().load(**kwargs)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/core/dataset.py", line 662, in load
    v.load()
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/core/variable.py", line 439, in load
    self._data = np.asarray(self._data)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/numpy/core/_asarray.py", line 83, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/core/indexing.py", line 685, in __array__
    self._ensure_cached()
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/core/indexing.py", line 682, in _ensure_cached
    self.array = NumpyIndexingAdapter(np.asarray(self.array))
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/numpy/core/_asarray.py", line 83, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/core/indexing.py", line 655, in __array__
    return np.asarray(self.array, dtype=dtype)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/numpy/core/_asarray.py", line 83, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/core/indexing.py", line 560, in __array__
    return np.asarray(array[self.key], dtype=None)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/numpy/core/_asarray.py", line 83, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/coding/variables.py", line 70, in __array__
    return self.func(self.array)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/xarray/coding/variables.py", line 220, in _scale_offset_decoding
    data *= scale_factor
ValueError: non-broadcastable output operand with shape () doesn't match the broadcast shape (1,)

What you expected to happen:

I expect the data access to work similarly as when opening with other engines.

Minimal Complete Verifiable Example:

import xarray
fn = "/data/gholl/cache/fogtools/abi/2017/03/14/20/06/7/OR_ABI-L1b-RadF-M3C07_G16_s20170732006100_e20170732016478_c20170732016514.nc"
with xarray.open_dataset(fn, engine="h5netcdf") as ds:
    ds["Rad"][400, 300].load()

Anything else we need to know?:

An earlier version of this issue, and some comments, refer to fsspec or working on open files, but that proved to have nothing to do with the problem.

Environment:

I've confirmed this issue installing xarray from latest master, which means xarray 0.16.2.dev11+gf821fe20 at the time of writing,

Output of xr.show_versions()

INSTALLED VERSIONS

commit: None
python: 3.8.5 | packaged by conda-forge | (default, Sep 24 2020, 16:55:52)
[GCC 7.5.0]
python-bits: 64
OS: Linux
OS-release: 4.12.14-lp150.12.82-default
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_GB.UTF-8
LOCALE: en_GB.UTF-8
libhdf5: 1.10.6
libnetcdf: 4.7.4

xarray: 0.16.2.dev11+gf821fe20
pandas: 1.1.2
numpy: 1.19.1
scipy: 1.5.2
netCDF4: 1.5.4
pydap: None
h5netcdf: 0.8.1
h5py: 2.10.0
Nio: None
zarr: 2.4.0
cftime: 1.2.1
nc_time_axis: None
PseudoNetCDF: None
rasterio: 1.1.6
cfgrib: None
iris: None
bottleneck: None
dask: 2.27.0
distributed: 2.27.0
matplotlib: 3.3.2
cartopy: None
seaborn: None
numbagg: None
pint: None
setuptools: 49.6.0.post20200917
pip: 20.2.3
conda: None
pytest: 6.0.2
IPython: None
sphinx: None

[djhoese]

Just tested this with decode_cf=False and the rest of the loading process seems fine (note: I used engine='h5netcdf' too).

[gerritholl]

I just tested this with some more combinations:

• decode_cf=True, mask_and_scale=False, everything seems fine.
• decode_cf=False, mask_and_scale=True, everything seems fine.
• decode_cf=True, mask_and_scale=True results in the ValueError and associated traceback.

[gerritholl]

Probably related: when reading an open file through a file system instance, the _FillValue, scale_factor, and add_offset are arrays of length one. When opening by passing a filename, those are all scalar (as expected):

import xarray
from fsspec.implementations.local import LocalFileSystem
fn = "/data/gholl/cache/fogtools/abi/2017/03/14/20/06/7/OR_ABI-L1b-RadF-M3C07_G16_s20170732006100_e20170732016478_c20170732016514.nc"
ds1 = xarray.open_dataset(fn, decode_cf=True, mask_and_scale=False)
print(ds1["esun"].attrs["_FillValue"])
print(ds1["Rad"].attrs["scale_factor"])
with LocalFileSystem().open(fn) as of:
    ds2 = xarray.open_dataset(of, decode_cf=True, mask_and_scale=False)
    print(ds2["esun"].attrs["_FillValue"])
    print(ds2["Rad"].attrs["scale_factor"])

Result:

-999.0
0.001564351
[-999.]
[0.00156435]

I strongly suspect that this is what causes the ValueError, and in any case it also causes downstream problems even if opening succeeds as per the previous comment.

[gerritholl]

Some further digging shows it's due to differences between the h5netcdf and netcdf4 backends:

import xarray
fn = "/data/gholl/cache/fogtools/abi/2017/03/14/20/06/7/OR_ABI-L1b-RadF-M3C07_G16_s20170732006100_e20170732016478_c20170732016514.nc"
with xarray.open_dataset(fn, decode_cf=False, mask_and_scale=False, engine="netcdf4") as ds:
    print(ds["esun"].attrs["_FillValue"])
    print(ds["Rad"].attrs["scale_factor"])
with xarray.open_dataset(fn, decode_cf=False, mask_and_scale=False, engine="h5netcdf") as ds:
    print(ds["esun"].attrs["_FillValue"])
    print(ds["Rad"].attrs["scale_factor"])

Results in:

-999.0
0.001564351
[-999.]
[0.00156435]

[gerritholl]

However, a simple `xarray.open_dataset(fn, engine="h5netcdf") still fails with ValueError only if passed an open file, so there appear to be still other differences apart from the dimensionality of the variable attributes depending on backend.

[gerritholl]

My last comment was inaccurate. Although the open succeeds, the non-scalar scale factor does trigger failure upon accessing data (due to lazy loading) even without any open file:

import xarray
fn = "OR_ABI-L1b-RadF-M3C07_G16_s20170732006100_e20170732016478_c20170732016514.nc"
with xarray.open_dataset(fn, engine="h5netcdf") as ds:
    print(ds["Rad"][400, 300])

The data file is publicly available at:

s3://noaa-goes16/ABI-L1b-RadF/2017/073/20/OR_ABI-L1b-RadF-M3C07_G16_s20170732006100_e20170732016478_c20170732016514.nc

[gerritholl]

Interestingly, the problem is prevented if one adds

ds.load()

before the print statement.

[gerritholl]

The ds.load() prevents the traceback because it means the entire n-d data variable is multiplied with the 1-d scale factor. Similarly, requesting a slice (ds["Rad"][400:402, 300:302]) also prevents a traceback. The traceback occurs if a single value is requested, because then Python will complain about multiplying a scalar with a 1-d array. I'm not entirely sure why, but would be a numpy issue:

In [7]: a = np.array(0)

In [8]: b = np.array([0])

In [9]: a * b
Out[9]: array([0])

In [10]: a *= b
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-10-0d04f348f081> in <module>
----> 1 a *= b

ValueError: non-broadcastable output operand with shape () doesn't match the broadcast shape (1,)