xarray: Problem opening unstructured grid ocean forecasts with 4D vertical coordinates

I’ve seen this issue come up a few more times recently, so I wanted to ask: in lieu of a “full fix” (a la https://github.com/pydata/xarray/pull/2405, with deep data model changes holding up the PR), would there be support for a partial coordinate-renaming-based fix? I’d imagine something like the following:

To read in netCDF like the following,

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    Conventions: CDM-Extended-CF
    history: Written by CFPointWriter
    title: Extract Points data from Grid file /data/ldm/pub/native/grid/NCEP/RR/CONUS_13km/RR_CONUS_13km_20210219_1400.grib2.ncx3#LambertConformal_337X451-40p01N-98p14W
    featureType: timeSeriesProfile
    time_coverage_start: 2021-02-19T16:00:00Z
    time_coverage_end: 2021-02-19T16:00:00Z
    geospatial_lat_min: 42.0295
    geospatial_lat_max: 42.0305
    geospatial_lon_min: -93.6405
    geospatial_lon_max: -93.6395
    dimensions(sizes): profile(1), station(1), isobaric(37), station_name_strlen(10), station_description_strlen(34)
    variables(dimensions): float32 isobaric(station, profile, isobaric), float32 u-component_of_wind_isobaric(station, profile, isobaric), float32 v-component_of_wind_isobaric(station, profile, isobaric), float32 Temperature_isobaric(station, profile, isobaric), float32 Relative_humidity_isobaric(station, profile, isobaric), |S1 station_name(station, station_name_strlen), |S1 station_description(station, station_description_strlen), float64 latitude(station), float64 longitude(station), float64 time(station, profile)
    groups: 

(note the problematic isobaric(station, profile, isobaric)), one could specify a kwarg to xr.open_dataset to rename it,

ds = xr.open_dataset("my_problematic_file.nc", rename_vars={'isobaric': 'isobaric_coord'})

thus giving

<xarray.Dataset>
Dimensions:                       (isobaric: 37, profile: 1, station: 1)
Dimensions without coordinates: isobaric, profile, station
Data variables:
    u-component_of_wind_isobaric  (station, profile, isobaric) float32 ...
    v-component_of_wind_isobaric  (station, profile, isobaric) float32 ...
    Temperature_isobaric          (station, profile, isobaric) float32 ...
    Relative_humidity_isobaric    (station, profile, isobaric) float32 ...
    station_name                  (station) |S10 ...
    station_description           (station) |S34 ...
    latitude                      (station) float64 ...
    longitude                     (station) float64 ...
    time                          (station, profile) datetime64[ns] ...
    isobaric_coord                (station, profile, isobaric) float32 ...
Attributes:
    Conventions:          CDM-Extended-CF
    history:              Written by CFPointWriter
    title:                Extract Points data from Grid file /data/ldm/pub/na...
    featureType:          timeSeriesProfile
    time_coverage_start:  2021-02-19T16:00:00Z
    time_coverage_end:    2021-02-19T16:00:00Z
    geospatial_lat_min:   42.0295
    geospatial_lat_max:   42.0305
    geospatial_lon_min:   -93.6405
    geospatial_lon_max:   -93.6395

This is the second follow-up item in https://github.com/pydata/xarray/issues/6293

I think we could definitely experiment with relaxing this constraint now, although ideally we would continue to check off auditing all of the methods in that long list first.

I had to go around this issue and not use xarray but pandas instead or plain netdcf4

nc = netCDF4.Dataset(input_file) h = nc.variables[vname] times = nc.variables[‘time’] jd = netCDF4.num2date(times[:],times.units) hs = pd.Series(h[:,station],index=jd)

I would love to know if there is a way to do it over xarray since it is so nice to use. Best regards

It is not ideal but you can workaround that by dropping the siglay variable.

ds = xr.open_dataset(url, drop_variables='siglay')

Hi

For now, I found a workaround loading and renaming the problematic coordinates with netCDF4.Dataset(). Soon I will post this and other solutions for this model output in iuryt/FVCOMpy.

For now, you could try:

import xarray as xr
from netCDF4 import Dataset

# define year and month to be read
year = 2019
month = 5
# we could use this to run a loop through the years/months we need


# list problematic coordinates
drop_variables = ['siglay','siglev']

# base url for openDAP server
url = "".join(["http://www.smast.umassd.edu:8080/thredds/dodsC/models/fvcom/",
        f"NECOFS/Archive/NECOFS_GOM/{year}/gom4_{year}{month:02d}.nc?"])

# lazy load of the data
ds = xr.open_dataset(url,drop_variables=drop_variables,decode_times=False)

# load data with netCDF4
nc = Dataset(url)
# load the problematic coordinates
coords = {name:nc[name] for name in drop_variables}
# function to extract ncattrs from `Dataset()`
get_attrs = lambda name: {attr:coords[name].getncattr(attr) for attr in coords[name].ncattrs()}
# function to convert from `Dataset()` to `xr.DataArray()`
nc2xr = lambda name: xr.DataArray(coords[name],attrs=get_attrs(name),name=f'{name}_coord',dims=(f'{name}','node'))
# merge `xr.DataArray()` objects
coords = xr.merge([nc2xr(name) for name in coords.keys()])
# reassign to the main `xr.Dataset()`
ds = ds.assign_coords(coords)

Leaving it here in case someone fall into the same problem.