Is there anyways to use HDF data using xarray?

[krishnaap]

I am wondering is there any way to use HDF files in xarray? For example, MOD06 files of the MODIS cloud dataset consist of multiple satellite swaths over a region given in multiple files, Is there anyway to open it using xarray and handle it as a xr.dataset?.

[dcherian]

There is no "open hdf5" but you can do it manually. Use h5py to open the file, extract the arrays you want, and pass that to xr.DataArray(...)

[krishnaap]

I tried it using multiple swaths (say 10 swaths over a wide region) to combine into a single file by saving into netCDF with xarray

import os
import glob                                                                 
import numpy as np
import xarray as xr

FILE_NAME = '/media/krishna/Linux/Dataset/MOD_APR_subsetted/HDF/MOD06_L2.A2019091.1555.061.2019092015305.hdf'
DATAFIELD_NAME = 'Cloud_Top_Pressure'

from pyhdf.SD import SD, SDC
#%%
i = 0

for file in list(glob.glob('MOD06_L2.A2019091*.hdf')):
    reader = open(file)
    hdf = SD(file, SDC.READ)
    # Read dataset.
    data2D = hdf.select(DATAFIELD_NAME)
    data = data2D[:,:].astype(np.double)
    # Read geolocation dataset.
    lat = hdf.select('Latitude')
    latitude = lat[:,:]
    lon = hdf.select('Longitude')
    longitude = lon[:,:]
    # Retrieve attributes.
    attrs = data2D.attributes(full=1)
    aoa=attrs["add_offset"]
    add_offset = aoa[0]
    fva=attrs["_FillValue"]
    _FillValue = fva[0]
    sfa=attrs["scale_factor"]
    scale_factor = sfa[0]        
    ua=attrs["units"]
    units = ua[0]
    data[data == _FillValue] = np.nan
    data = (data - add_offset) * scale_factor 
    datam = np.ma.masked_array(data, np.isnan(data))
    if i == 0 :
        data_m = datam
        latitude_m = latitude
        longitude_m = longitude
    else:
        data_m = np.vstack([data_m, datam])
        print(data_m)

        latitude_m = np.vstack([latitude_m, latitude])
        longitude_m = np.vstack([longitude_m, longitude])
    i = i + 1
latitude = latitude_m
longitude = longitude_m
ctp = data_m
# ds = xr.Dataset(data_m,longitude_m, latitude_m)

#%%
df = xr.Dataset(
        data_vars={'ctp':    (('latitude', 
'longitude' ), ctp)},
        coords={'longitude': longitude,
            'latitude': latitude,
        })
df.ctp.attrs['long_name'] = 'Cloud Top Pressure'


df.longitude.attrs['standard_name'] = 'longitude'
df.longitude.attrs['long_name'] = 'longitude'
df.longitude.attrs['units'] = 'degrees_east'
df.longitude.attrs['axis'] = 'X'

df.latitude.attrs['standard_name'] = 'latitude'
df.latitude.attrs['long_name'] = 'latitude'
df.latitude.attrs['units'] = 'degrees_north'
df.latitude.attrs['axis'] = 'Y'


df.to_netcdf('MOD_2019-04-01.nc')

But it's showing an error

MissingDimensionsError: cannot set variable 'longitude' with 2-dimensional data without explicit dimension names. Pass a tuple of (dims, data) instead.

Traceback (most recent call last):

File "/media/krishna/Linux/Dataset/MOD_APR_subsetted/HDF/MODIS_day_saving.py", line 97, in
'latitude': latitude,

File "/home/krishna/anaconda3/envs/unified/lib/python3.6/site-packages/xarray/core/dataset.py", line 742, in init
data_vars, coords, compat="broadcast_equals"

File "/home/krishna/anaconda3/envs/unified/lib/python3.6/site-packages/xarray/core/merge.py", line 476, in merge_data_and_coords
indexes = dict(_extract_indexes_from_coords(coords))

File "/home/krishna/anaconda3/envs/unified/lib/python3.6/site-packages/xarray/core/merge.py", line 485, in _extract_indexes_from_coords
variable = as_variable(variable, name=name)

File "/home/krishna/anaconda3/envs/unified/lib/python3.6/site-packages/xarray/core/variable.py", line 155, in as_variable
"(dims, data) instead." % (name, data.ndim)

MissingDimensionsError: cannot set variable 'longitude' with 2-dimensional data without explicit dimension names. Pass a tuple of (dims, data) instead.

[JessicaS11]

Depending on how many variables you want, you can also use one of the netcdf engines to read a group-path/variable at a time into DataSets, which you'll then have to manually manipulate a bunch to get coordinates and dimensions right and merge into a single DataSet.

[krishnaap]

I tried to open it using the netCDF4 engine too, but it is showing an error for multiple files.

ValueError: arguments without labels along dimension 'Cell_Along_Swath_1km:mod06' cannot be aligned because they have different dimension sizes: {2040, 2030}

MODIS swaths come in two sizes in L2 products, especially for aerosol and cloud products. Is there anyway to incorporate this into the code?

import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy as ccrs
import numpy as np

fname='/Dataset/MOD06_L2/MOD06_L2/2019/091/MOD06_L2.A2019091*.hdf'        

ds=xr.open_mfdataset(fname,engine='netcdf4',concat_dim='Cell_Across_Swath_hkm = 2708',combine='nested')
varlist=list(ds.variables)

#%%
ctp = ds.Cloud_Top_Pressure
lat = ds.Latitude
lon = ds.Longitude
ctp

#%%
plt.figure(figsize=(14, 6))
ax = plt.axes(projection=ccrs.crs.PlateCarree())
ax.set_title('Cloud Top Pressure (mb)')
ax.set_global()
ctp[0].plot.pcolormesh(ax=ax, transform=ccrs.crs.PlateCarree(),cmap=plt.cm.jet,label="Cloud Top Pressure over North East")#, add_colorbar=False)
ax.coastlines()
ax.set_xlim([60, 100]);
ax.set_ylim([0, 40]);
ax.set_yticks(np.arange(0, 40, 10))
ax.set_xticks(np.arange(60, 100, 10))

[krishnaap]

@JessicaS11 is there any solution for this? I am kind of stuck at this point

[rwegener2]

Hi there,

I'm coming back to this topic because I had a similar question about the recommended way to open a single hdf5 group in xarray.

On the Reading and Writing page of docs I noticed that HDF5 is primarily mentioned in the chunking section, and once in comparison to netCDF4. Using a sort of "guess and check" method I figured out that you can open a particular HDF5 group using the same syntax described in the netCDF4 groups section:

xr.open_dataset('path/to/file', group='/my/group')

Is that the recommended way to open a single HDF5 group in xarray? If so, would it make sense to add an "HDF5" section to the Reading and Writing docs page? I would be happy to submit a draft PR for review. Main points would be something like:

• HDF5 files can be opened using the h5netcdf engine
• xarray cannot interrogate an HDF5 file for group names, but if you know the name of the group you want the group can be specified with the group argument.
• others?

[jhamman]

I'm not surprised to hear that reading large and/or deeply nested hdf5 files in datatree is slow. We simply haven't done anything to optimize the performance in this part of the library yet. @jthielen got a start on a native datatree backend but I think that effort has stalled. Perhaps it could be resurrected?

[TomNicholas]

It was this PR: #7437

[keewis]

Perhaps it could be resurrected?

definitely, and I'm currently working on something that would be easier if the backend system would exist already, so I'm interested in helping out.

[rwegener2]

Very nice that there is already a PR going!

For my understanding, @jhamman, why would adding datatree as a native xarray backend speed up reads?

[jthielen]

Perhaps it could be resurrected?

definitely, and I'm currently working on something that would be easier if the backend system would exist already, so I'm interested in helping out.

Sorry for the delay, but I'd be glad to work on this PR again! I'd appreciate a review of the general approach taken to implement it, before moving on (to things like tests and de-duplicating some code).