How do I open datasets with nested groups using xr.open_mfdatasets?

[mccabete]

Description

I am working with precipitation data. I was delighted that some precip products, like the daily IMERG '10.5067/GPM/IMERGDF/DAY/07' work out-of-the-box using the provided cloud-access documentation example.

Other products are more complicated. For example, I figured out that I needed to include more parameters in xr.open_mfdataset() for half-hourly ('10.5067/GPM/IMERG/3B-HH/07'), because it sotred within the group "grid".

ds = xr.open_mfdataset(fn, group="Grid", concat_dim="time", combine="nested", coords="all", data_vars="all")

But I am still struggling with how to open files with a more complicated nesting structure, like the GPROF data ('10.5067/GPM/GMI/GPM/GPROF/2A/07'), which is nested at multiple levels. Is there a standard protocol?

[betolink]

For really nested datasets xarray has open_datatree(). This method works at the single file level, you can see the structure and dimensions of the nested datasets/data-arrays but as far as I know it doesn't have a "open many of this" implementation like open_mfdataset. Perhaps we can use open_datatree() in the preprocess step of open_mfdataset. I can take a look and will report back.

Do you guys have tutorials on this dataset @battistowx?

[jhkennedy]

I asked why there was no open_mfdatatree recently on Earthmover's slack and these were the important responses:

I wouldn't want to build it because it'd take a ton of work to get remotely acceptable performance, and still wouldn't be as efficient as using a open_datatree on a virtual zarr store

in theory we could, but IMO the functions' signature would likely either be hugely complicated or too opinionated to be useful

It would be a god function. Kinda the whole point of VirtualiZarr is that there isn't just one god function, as you have xarray's primitives of open/concat/merge to build up what you actually need

So the recommendation is to use VirtualiZarr to build a virtual "view" into the datasets with everything you're interested in.

[mccabete]

This is great feedback! The next question is "how does someone use VirtualliZarr to 'view' these nested datasets?". It looked like the earthaccess community already has a similar example here, where earthacess.visualize() is used for a dataset that has multiple top-level groups/ a clear dimention that the datasets can be aggrigated over. But I can't immediately see how one would use it for nested groups, or for a dataset where those nested groups hold the dimension I'd like to concatenate over (time). Trying out variations of that code throws what looks like different structural errors:

import earthaccess
import numpy as np
import xarray as xr

auth = earthaccess.login()

results = earthaccess.search_data(
    doi = '10.5067/GPM/GMI/GPM/GPROF/2A/07', 
    cloud_hosted=True,
    temporal=("2023-07-28T00:00:00", "2023-07-31T12:00:00"),
    point = (-115.22265, 35.25335)
)

open_options = {
    "access": "indirect", 
    "load": True,  
    "concat_dim": "time",  
    "data_vars": "minimal",  
    "coords": "minimal",  
    "compat": "override",  
    "combine_attrs": "override", 
}

result_root = earthaccess.virtualize(granules=results, **open_options) ## Fails with 'ValueError: Cannot specify both data_vars='minimal' and coords='minimal' when concatenating over a new dimension'

result_s1 = earthaccess.virtualize(
    granules=results, group="S1", **open_options
) # Fails with 'KeyError: 'name'

result_header = earthaccess.virtualize(
    granules=results, group="GprofDHeadr", **open_options
) # Fails with 'KeyError: 'name'

For reference, the structure of TEMPO data the example works with is

geolocation
geolocation/relative_azimuth_angle
geolocation/solar_zenith_angle
geolocation/viewing_zenith_angle
latitude
longitude
product
product/main_data_quality_flag
product/vertical_column_stratosphere
product/vertical_column_troposphere
product/vertical_column_troposphere_uncertainty
qa_statistics
qa_statistics/max_vertical_column_stratosphere_sample
qa_statistics/max_vertical_column_total_sample
qa_statistics/max_vertical_column_troposphere_sample
qa_statistics/max_vertical_column_troposphere_uncertainty_sample
qa_statistics/min_vertical_column_stratosphere_sample
qa_statistics/min_vertical_column_total_sample
qa_statistics/min_vertical_column_troposphere_sample
qa_statistics/min_vertical_column_troposphere_uncertainty_sample
qa_statistics/num_vertical_column_stratosphere_samples
qa_statistics/num_vertical_column_total_samples
qa_statistics/num_vertical_column_troposphere_samples
qa_statistics/num_vertical_column_troposphere_uncertainty_samples
support_data
support_data/albedo
support_data/amf_cloud_fraction
support_data/amf_cloud_pressure
support_data/amf_stratosphere
support_data/amf_total
support_data/amf_troposphere
support_data/eff_cloud_fraction
support_data/fitted_slant_column
support_data/fitted_slant_column_uncertainty
support_data/snow_ice_fraction
support_data/surface_pressure
support_data/terrain_height
support_data/tropopause_pressure
support_data/vertical_column_total
support_data/vertical_column_total_uncertainty
time
weight

```



But the structure for the GPROF data looks like this: 




```

GprofDHeadr
GprofDHeadr/clusterProfiles
GprofDHeadr/hgtTopLayer
GprofDHeadr/speciesDescription
GprofDHeadr/temperatureDescriptions
S1
S1/L1CqualityFlag
S1/Latitude
S1/Longitude
S1/SCstatus
S1/SCstatus/FractionalGranuleNumber
S1/SCstatus/SCaltitude
S1/SCstatus/SClatitude
S1/SCstatus/SClongitude
S1/SCstatus/SCorientation
S1/ScanTime
S1/ScanTime/DayOfMonth
S1/ScanTime/DayOfYear
S1/ScanTime/Hour
S1/ScanTime/MilliSecond
S1/ScanTime/Minute
S1/ScanTime/Month
S1/ScanTime/Second
S1/ScanTime/SecondOfDay
S1/ScanTime/Year
S1/airmassLiftIndex
S1/cloudWaterPath
S1/convectivePrecipitation
S1/frozenPrecipitation
S1/iceWaterPath
S1/mostLikelyPrecipitation
S1/pixelStatus
S1/precip1stTertial
S1/precip2ndTertial
S1/precipitationYesNoFlag
S1/probabilityOfPrecip
S1/profileNumber
S1/profileScale
S1/profileTemp2mIndex
S1/qualityFlag
S1/rainWaterPath
S1/sunGlintAngle
S1/sunLocalTime
S1/surfacePrecipitation
S1/surfaceTypeIndex
S1/temp2mIndex
S1/totalColumnWaterVaporIndex

[battistowx]

For really nested datasets xarray has open_datatree(). This method works at the single file level, you can see the structure and dimensions of the nested datasets/data-arrays but as far as I know it doesn't have a "open many of this" implementation like open_mfdataset. Perhaps we can use open_datatree() in the preprocess step of open_mfdataset. I can take a look and will report back.

Do you guys have tutorials on this dataset @battistowx?

I don't have this published anywhere else, but this is a Jupyterlite notebook showing how to use open_datatree on an example GPM collection: https://cbattisto-demos.github.io/DataTree-Demo/lab/index.html.

We are actually still in the process of testing if GPM IMERG (Level 3) can be virtualized, and there are many caveats with attempting to virtualize Level 2 collections like GPM_GMI in the example code from @mccabete, which we have not tested yet. I would recommend trying the datatree solution before virtualization.

[itcarroll]

I just want to make sure it has been noted that the group argument in open_mfdataset accepts a "/" seperated path to nested groups. You can do xr.open_mfdataset(..., group="S1/SCstatus") for instance. Dealing with coords in other groups remains a problem.