HDF error when trying to write Dataset read with rasterio to NetCDF

[loicdtx]

I'm getting an HDF error when trying to write a Dataset read from GeoTiff (rasterio backend) to NetCDF. See reproducible example below:

import urllib.request
import tempfile
import os

import xarray as xr

path = tempfile.gettempdir()
url = 'https://earthexplorer.usgs.gov/browse/gisready/landsat_8/LC08_L1TP_026047_20180110_20180119_01_T1.zip'
filename = os.path.join(path, url.split('/')[-1])
nc_name = os.path.join(path, 'landsat_rgb.nc')

# Download file if not exist (11 Mb)
if not os.path.isfile(filename):
    urllib.request.urlretrieve(url, filename)

# Read rgb file using rasterio backend
rgb_name = '/'.join(['/vsizip', filename,
                     os.path.basename(filename).split('.')[-2] + '.tif'])
ds = xr.open_rasterio(rgb_name)
ds = ds.to_dataset('band').rename({1:'blue', 2:'green', 3:'red'})
print(ds)

# <xarray.Dataset>
# Dimensions:  (x: 7611, y: 7761)
# Coordinates:
#   * y        (y) float64 2.193e+06 2.193e+06 2.193e+06 ... 1.961e+06 1.960e+06
#   * x        (x) float64 3.732e+05 3.732e+05 3.733e+05 ... 6.015e+05 6.015e+05
# Data variables:
#     blue     (y, x) uint8 ...
#     red      (y, x) uint8 ...
#     green    (y, x) uint8 ...
# Attributes:
#     transform:   (30.0, 0.0, 373185.0, 0.0, -30.0, 2193315.0)
#     crs:         +init=epsg:32614
#     res:         (30.0, 30.0)
#     is_tiled:    1
#     nodatavals:  (nan, nan, nan)


# Write to netcdf
ds.to_netcdf(nc_name)

Output of xr.show_versions()

python -c "import xarray as xr; xr.show_versions()"

INSTALLED VERSIONS

commit: None
python: 3.6.5.final.0
python-bits: 64
OS: Linux
OS-release: 4.15.0-36-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: en_US.UTF-8

xarray: 0.10.9
pandas: 0.23.4
numpy: 1.15.3
scipy: None
netCDF4: 1.4.2
h5netcdf: None
h5py: None
Nio: None
zarr: 2.2.0
cftime: 1.0.2.1
PseudonetCDF: None
rasterio: 1.0.9
iris: None
bottleneck: None
cyordereddict: None
dask: 0.20.0
distributed: None
matplotlib: None
cartopy: None
seaborn: None
setuptools: 40.5.0
pip: 18.1
conda: None
pytest: None
IPython: 7.1.1
sphinx: None

[jhamman]

@loicdtx - thanks for the report. I just tried this using xarray master and didn't get an error. So, you could try that and see if that works. Since our last release, we've had a fairly significant refactor of the IO backends in xarray.

[loicdtx]

Hi @jhamman, just tried what you suggest. Apparently it's a scipy vs netCDF4 thing; it works with scipy but not with netCDF4 on both master HEAD and the latest stable release.

[jhamman]

Interesting. I'm not seeing any difference when using scipy or netcdf4 or h5netcdf.

(by the way, thank you for the reproducible example)

[loicdtx]

yes, that's interesting... are you using a different OS? Let me know if there's something else I can help with.

[shoyer]

Exactly what error message do you see?

[loicdtx]

@shoyer, below the full traceback; I installed everything with pip

mktmpenv -p python3
pip install xarray numpy netcdf4
pip install rasterio

Traceback (most recent call last):
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/api.py", line 724, in to_netcdf
    unlimited_dims=unlimited_dims, compute=compute)
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/core/dataset.py", line 1179, in dump_to_store
    unlimited_dims=unlimited_dims)
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/common.py", line 374, in store
    unlimited_dims=unlimited_dims)
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/netCDF4_.py", line 406, in set_variables
    super(NetCDF4DataStore, self).set_variables(*args, **kwargs)
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/common.py", line 413, in set_variables
    self.writer.add(source, target)
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/common.py", line 272, in add
    target[...] = source
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/netCDF4_.py", line 48, in __setitem__
    data[key] = value
  File "netCDF4/_netCDF4.pyx", line 4648, in netCDF4._netCDF4.Variable.__setitem__
  File "netCDF4/_netCDF4.pyx", line 4913, in netCDF4._netCDF4.Variable._put
  File "netCDF4/_netCDF4.pyx", line 1754, in netCDF4._netCDF4._ensure_nc_success
RuntimeError: NetCDF: HDF error

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "rasterio_reprex.py", line 41, in <module>
    ds.to_netcdf(nc_name)
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/core/dataset.py", line 1254, in to_netcdf
    compute=compute)
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/api.py", line 729, in to_netcdf
    store.close()
  File "/home/loic/.virtualenvs/tmp-3caa6b25124e2f31/lib/python3.6/site-packages/xarray/backends/netCDF4_.py", line 474, in close
    ds.close()
  File "netCDF4/_netCDF4.pyx", line 2276, in netCDF4._netCDF4.Dataset.close
  File "netCDF4/_netCDF4.pyx", line 2260, in netCDF4._netCDF4.Dataset._close
  File "netCDF4/_netCDF4.pyx", line 1754, in netCDF4._netCDF4._ensure_nc_success
RuntimeError: NetCDF: HDF error

[jhamman]

Is this the same error you get with xarray/master?

[fmaussion]

I can reproduce this on my machine - also linux. This is going to be hard to track down though.

[loicdtx]

@jhamman, yes, same error message when installing from master

[ghost]

I have similar problem, when importing rasterio in the same script (not even using it for anything). This fails with HDF error:

import xarray as xa
import numpy as np
#import netCDF4
import rasterio

ds = xa.Dataset()
ds['z'] = (('y', 'x'), np.zeros((100, 100), np.float32))
print(ds)
ds.to_netcdf('test.nc')
ds.close()

with xa.open_dataset('test.nc') as ds:
    print(ds)

If I import netCDF4 before rasterio it works fine (uncomment line 3). This is probably an issue with rasterio somehow.

I installed everything with pip:

$ pip install Cython
$ pip install netCDF4 xarray rasterio numpy

From pip freeze:

affine==2.2.1
attrs==18.2.0
cftime==1.0.3
Click==7.0
click-plugins==1.0.4
cligj==0.5.0
Cython==0.29.1
netCDF4==1.4.2
numpy==1.15.4
pandas==0.23.4
pyparsing==2.3.0
python-dateutil==2.7.5
pytz==2018.7
rasterio==1.0.11
six==1.11.0
snuggs==1.4.2
xarray==0.11.0

[ghost]

It seems that this is not a problem with xarray but only with rasterio and netCDF4. Also this fails:

import rasterio
import netCDF4

with netCDF4.Dataset('test.nc', mode='w') as ds:
    ds.createDimension('x')
    ds.createVariable('foo', float, dimensions=('x'))
    print(ds)

Commenting out import rasterio removes the HDF error. I’ll report this to rasterio.

[shoyer]

This looks like a binary incompatibility issue with wheels for rasterio and netCDF4 on PyPI.

Good alternatives (for now) would be building from source or using conda-forge.

[ChristianF88]

Hi there,

I just wanted to let you know, that I do get the same error, when working with the following script. It does not use xarray or rasterio. So those will most likely not be the problem.

class Gif2NetCDF():
    def __init__(self, netcdf_name, netcdf_destfolder, gif_folder,
                 gif_filter_pattern=None, detailled_conversion=True):
        """

        Parameters
        ----------
        netcdf_name : str
            specifying the name of the netcdf file to write.

        netcdf_destfolder : str
            specifying the folder where the netcdf file is supposed to be created

        giffolder : str
            specifying the folder that contains the gif files, that are supposed to be written as an netcdf

        gif_filter_pattern : str
            specifying a re pattern to filter the all files in directory so you end up with the gif files you want.
            The here specified string will be passed to re.match which will be checked for each file in the provided giffolder


        Examples
        --------

        ## defining variables
        import netCDF4
        netcdf_name = "2016-07.nc"
        netcdf_destfolder = "./example/radar"
        giffolder = "./example/radar/2016-07/"

        ## create Instance
        cdf = Gif2NetCDF(netcdf_name, netcdf_destfolder, giffolder,gif_filter_pattern)
        ## write all (filtered) gifs in folder to netcdf file
        cdf.writeCDF()
        """

        ## creating global vars

        self.netcdf_name = netcdf_name
        self.netcdf_destfolder = netcdf_destfolder
        self.giffolder = gif_folder
        self.refilterpattern = gif_filter_pattern
        self.detailled_conversion = detailled_conversion

        self.netcdfFP = os.path.join(netcdf_destfolder, netcdf_name)

        # preparing the coordinate vectors
        self._lat_range = np.arange(479500, -160500, step=-1000)  # north
        self._long_range = np.arange(255500, 965500, step=1000)  # east

        # preparing time origin
        self.time_origin = dt.datetime.strptime("1970-01-01 00:00:00", "%Y-%m-%d %H:%M:%S")
        
        self.raincodes = np.array([....]) # this array is quite large, so I left it out... Which unfortunately means the code does not run... If anybody needs it please let me know

        return

    def list_gifs(self):
        self._gifs = os.listdir(self.giffolder)
        return self._gifs

    def filter_gifs(self):
        self._gifs = [file for file in self._gifs if re.match(self.refilterpattern, file)]
        return self._gifs

    def addDimensions(self):

        # adds dimensions to empty netcdf file
        self.latD = self.netcdfFile.createDimension('lat', self._lat_range.shape[0])  # north-south
        self.lonD = self.netcdfFile.createDimension('lon', self._long_range.shape[0])  # east-west
        self.timeD = self.netcdfFile.createDimension('time', None)

        return

    def addVariables(self):

        ## creating variables
        self.latV = self.netcdfFile.createVariable("chy", np.float32, ("lat",), complevel=9, zlib=True)  # north-south
        self.lonV = self.netcdfFile.createVariable("chx", np.float32, ("lon",), complevel=9, zlib=True)  # east-west
        self.timeV = self.netcdfFile.createVariable("time", np.float64, ("time",), complevel=9, zlib=True)

        self.rainV = self.netcdfFile.createVariable("rain", np.float32, ("time", "lat", "lon"), complevel=9, zlib=True,
                                                    fill_value=-100)

        ## adding units
        self.latV.units = "meters"
        self.lonV.units = "meters"

        self.timeV.units = "seconds"
        self.timeV.calender = "standard"

        self.rainV.units = "millimeter/hour"

        ## adding longname
        self.latV.long_name = "swiss northing CH1903"
        self.lonV.long_name = "swiss easting CH1903"
        self.timeV.long_name = "seconds since 1970-01-01 00:00:00"

        self.rainV.long_name = "precipitation intensity forecast"

        return

    def addDescription(self):

        self.netcdfFile.description = """..."""

        self.netcdfFile.history = """Created: {}""".format(dt.datetime.now().strftime("%Y-%m-%d %H:%M"))
        self.netcdfFile.source = '...'
        return

    def _write_static_dims(self):
        self.latV[:] = self._lat_range
        self.lonV[:] = self._long_range
        return

    def _write_time(self, file, datetime=None):

        if datetime is None:
            datestr = re.findall("\.([0-9]+)\.gif", file)[0]
            date = dt.datetime.strptime(datestr, "%Y%m%d%H%M")
        else:
            date = datetime

        seconds = (date - self.time_origin).total_seconds()
        current_size = self.timeV.size
        self.timeV[current_size] = seconds
        
        return

    def gif2array(self, file):
    
        xpix = 0
        ypix = 76
        n_pixel_x = 710 + xpix
        n_pixel_y = 640 + ypix

        gif = np.array(Image.open(file))[ypix:n_pixel_y, xpix:n_pixel_x].astype("float64")
        for idx, raincode in enumerate(self.raincodes):
            gif[gif == idx] = raincode[3]
        return gif

    def _write_rain(self, file):
        array = self.gif2array(os.path.join(self.giffolder, file))
        idx = self.rainV.shape[0] - 1
        self.rainV[idx, :, :] = array
        return

    def writeCDF(self):
        self.netcdfFile = Dataset(self.netcdfFP, 'w', format='NETCDF4_CLASSIC', )
        try:
            giflist = self.list_gifs()
            if self.refilterpattern is not None:
                fgiflist = self.filter_gifs()

            self.addDimensions()
            self.addVariables()
            self.addDescription()
            self._write_static_dims()

            for file in tqdm(self._gifs):
                self._write_time(file)
                self._write_rain(file)
                
        except Exception:
            self.netcdfFile.close()
            raise

        self.netcdfFile.close()

        return

Error

Traceback (most recent call last):
.
.
.
  File "C:\Users\foerstch\AppData\Local\Programs\Python\Python37\lib\site-packages\archiving\radar.py", line 358, in _write_rain
    idx = self.rainV.shape[0] - 1
  File "netCDF4\_netCDF4.pyx", line 4031, in netCDF4._netCDF4.Variable.shape.__get__
  File "netCDF4\_netCDF4.pyx", line 3369, in netCDF4._netCDF4.Dimension.__len__
  File "netCDF4\_netCDF4.pyx", line 1857, in netCDF4._netCDF4._ensure_nc_success
RuntimeError: NetCDF: HDF error

Does anybody have some advice on how to fix this?

Thanks a bunch!
Christian

[dcherian]

HDF error usually means a corrupt file. Does ncdump -h work on your file?

I've found that sometimes ncdump -h will succeed but there'll still be some corrupt data which will result in an error when you try to load variable values

[ChristianF88]

Well the error occurs while writing the file. So if the reason is, that the file is corrupt then necdf4 is corrupting it... right? The problem is rather, that it does not even finish writing the file.

FYI: The error does not occur at the same file every time.

Have a good weekend!! ;)

[alpha-beta-soup]

Similar error when using xr.open_rasterio, a workaround seems to be to change the order in which my datasets are opened. Example:

import xarray as xr
ds = xr.open_dataset('/data/someFile.nc') # netcdf
m = xr.open_rasterio('/data/otherFile.tif') # geotif
# Everything is happy

import xarray as xr
m = xr.open_rasterio('/data/otherFile.tif') # geotif
ds = xr.open_dataset('/data/someFile.nc') # netcdf
# Results in the following error

Traceback (most recent call last):
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/file_manager.py", line 186, in _acquire_with_cache_info
    file = self._cache[self._key]
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/lru_cache.py", line 42, in __getitem__
    value = self._cache[key]
KeyError: [<class 'netCDF4._netCDF4.Dataset'>, ('/data/someFile.nc',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False))]

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/api.py", line 420, in open_dataset
    filename_or_obj, group=group, lock=lock, **backend_kwargs)
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/netCDF4_.py", line 335, in open
    autoclose=autoclose)
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/netCDF4_.py", line 293, in __init__
    self.format = self.ds.data_model
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/netCDF4_.py", line 344, in ds
    return self._acquire()
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/netCDF4_.py", line 338, in _acquire
    with self._manager.acquire_context(needs_lock) as root:
  File "/usr/lib/python3.6/contextlib.py", line 81, in __enter__
    return next(self.gen)
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/file_manager.py", line 174, in acquire_context
    file, cached = self._acquire_with_cache_info(needs_lock)
  File "/usr/local/lib/python3.6/dist-packages/xarray/backends/file_manager.py", line 192, in _acquire_with_cache_info
    file = self._opener(*self._args, **kwargs)
  File "netCDF4/_netCDF4.pyx", line 2291, in netCDF4._netCDF4.Dataset.__init__
  File "netCDF4/_netCDF4.pyx", line 1855, in netCDF4._netCDF4._ensure_nc_success
OSError: [Errno -101] NetCDF: HDF error: b'/data/someFile.nc'

>>> xr.show_versions()

INSTALLED VERSIONS
------------------
commit: None
python: 3.6.6 (default, Sep 12 2018, 18:26:19) 
[GCC 8.0.1 20180414 (experimental) [trunk revision 259383]]
python-bits: 64
OS: Linux
OS-release: 4.15.0-58-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: C.UTF-8
LANG: C.UTF-8
LOCALE: en_US.UTF-8
libhdf5: 1.10.2
libnetcdf: 4.6.3

xarray: 0.12.3
pandas: 0.25.1
numpy: 1.13.3
scipy: 1.3.1
netCDF4: 1.5.1.2
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.0.3.4
nc_time_axis: None
PseudoNetCDF: None
rasterio: 1.0.26
cfgrib: None
iris: None
bottleneck: None
dask: 2.3.0
distributed: 2.3.2
matplotlib: 3.1.1
cartopy: 0.17.0
seaborn: None
numbagg: None
setuptools: 39.0.1
pip: 9.0.1
conda: None
pytest: None
IPython: None
sphinx: None

And I have the HDF5_USE_FILE_LOCKING environment variable set to FALSE.

[kmuehlbauer]

Just a heads up, I can't reproduce this with latest packages from conda-forge for neither for xr.open_rasterio nor xr.open_dataset with engine=rasterio. Maybe the issue got resolved upstream.