Accessing 2D cross-sections of a large 4D dataset

[velimskyj]

Hello,
I have a large 4-D dataset in NetCDF4 file, created by NetCDF4 library in Fortran. I use level 1 compression in NetCDF creation, and it looks like this:

~:ncinfo JLphi_run_011a.nc
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
dimensions(sizes): longitude(720), latitude(360), depth(75), time(3652)
variables(dimensions): float64 longitude(longitude), float64 latitude(latitude), float64 depth(depth), int32 time(time), float32 JL_phi(time,depth,latitude,longitude)
groups:

When I import the file into an xarray and plot 2-D lat/lon cross-sections for different depths and times, everything works fine:

drange= [1.0, 2.0, 5.0, 10.0, 20.0, 50.0 ]
trange= ['2008-03-21T12:00:00', '2008-06-21T12:00:00', '2008-09-23T12:00:00', '2008-12-21T12:00:00']

X= xr.open_dataarray(InFile)
for t in trange:
  for d in drange:
    X2d= X.sel(time=t,depth=d,method='nearest')
    X2d.plot(...)

But if I want to plot 2-D time/depth plots for fixed geographic coordinates, the execution gets extremely slow
and often fails with RuntimeError: NetCDF: HDF error:

lat_points = np.array([   35, -74, 17, 70, 35, -30])
lon_points = np.array([285.5, 306, 90, 40, 20, 200])
X= xr.open_dataarray(InFile)
X2d= X.sel(latitude=xr.DataArray(lat_points),longitude=xr.DataArray(lon_points),method='nearest')
for l in range(npoints):
  X2d.isel(dim_0=l).plot(...)

I have experimented with different engine settings in open_dataarray as well as with the .load() applied
to X or X2d, but it didn't help. My Fortran experience tells me that the array ordering in memory is the key,
but I do not know enough of xarray internals. I would welcome any advise. Thank you.

[dcherian]

The slowness is expected and fits with your array ordering intuition.

HDF Error is not Xarray's fault; I'd look in to hdf5 logs somehow to find the source of the error.