Error reading variables from NetCDF file with Xarray when there are multiple variable layers

Question

I get an error reading variables from NetCDF file with Xarray when there are multiple variable layers

Data:

The netCDF file that I have has multiple layers:

• /Geophysical_Data/precipitation_total_surface_flux
• /cell_lat
• /cell_row

Issue:

When I read this file with Xarray, the data one layer below (i.e., /Geophysical_Data/precipitation_total_surface_flux) won't be read as variables, but as attributes.

import xarray as xr
ds = xr.open_data(file_path)

This happens when I use a default 'netcdf4' engine. When I change the engine to 'rasterio' the data is successfully read, but it takes more time to load + I get another issue where Y coordinate is read wrong, so I want to avoid it.

Therefore, I am looking for ways to read data with Xarray + netcdf4 engine by specifying layers to read (here /Geophysical_Data/.).

Note that I downloaded this data via the NSIDC server and I subset it there before download.

Snipped:

With netcdf4 engine

<xarray.Dataset>
Dimensions:                 (y: 1822, x: 3856)
Coordinates:
  * x                       (x) float64 -180.0 -179.9 -179.8 ... 179.9 180.0
  * y                       (y) float64 85.0 84.9 84.81 ... -84.83 -84.92 -85.01
Data variables:
    cell_column             (y, x) float64 ...
    cell_lat                (y, x) float32 ...
    cell_lon                (y, x) float32 ...
    cell_row                (y, x) float64 ...
    projection_information  object ...
    time                    datetime64[ns] ...
Attributes:
    Source:        v17.11.1
    Institution:   NASA Global Modeling and Assimilation Office
    History:       File written by ldas2daac.x
    Comment:       HDF-5
    Filename:      /discover/nobackup/projects/gmao/smap/SMAP_L4/L4_SM/Vv7032...
    Title:         SMAP L4_SM Geophysical (GPH) Data Granule
    grid_rows:     1822
    grid_columns:  3856
    Conventions:   CF
    References:    see SMAP L4_SM Product Specification Documentation
    Contact:       http://gmao.gsfc.nasa.gov

With rasterio engine:

<xarray.Dataset>
Dimensions:                           (band: 1, x: 3856, y: 1822)
Coordinates:
  * band                              (band) int32 1
  * x                                 (x) float64 -180.0 -179.9 ... 179.9 180.0
  * y                                 (y) float64 85.0 84.9 ... -84.92 -85.01
    projection_information            int32 ...
Data variables:
    cell_column                       (band, y, x) float64 ...
    cell_lat                          (band, y, x) float32 ...
    cell_lon                          (band, y, x) float32 ...
    cell_row                          (band, y, x) float64 ...
    precipitation_total_surface_flux  (band, y, x) float32 ...
Attributes: (12/201)
    antennaRotationRate:                13
    description:                        The observation-corrected surface met...
    identifier:                         GEOS5COR
    edition:                            JPL CL#14-2285, JPL 400-1567
    publicationDate:                    2014-07-01
    title:                              SMAP Handbook
    ...                                 ...
    grid_rows:                          1822
    History:                            File written by ldas2daac.x
    Institution:                        NASA Global Modeling and Assimilation...
    References:                         see SMAP L4_SM Product Specification ...
    Source:                             v17.11.1
    Title:                              SMAP L4_SM Geophysical (GPH) Data Gra...

Answer 1

Duplicate to this question & answer. group argument in open_dataset() method will do the job.

To merge data from different layers, do the followings (referring to kmuehlbauer's answer on the Github issue).

# load root group with coordinates
ds_root = xr.open_dataset(os.path.join(input_path, fn), group="/", engine='netcdf4')
# load data from Geophysical_Data group
ds_precip = xr.open_dataset(os.path.join(input_path, fn), group="Geophysical_Data", engine='netcdf4')
# merge groups
ds = xr.merge([ds_root , ds_precip])
# plot
ds .precipitation_total_surface_flux.plot()