https://github.com/pydata/xarray/issues/2233#issuecomment-397602084 Would the new xarray index/coordinate internal refactoring now allow us to address this issue?

cc @kthyng

@benbovy Ah, I see you mean under "Relax all constraints related to “dimension (index) coordinates” in Xarray". Ok, thank you for clarifying that for me! (I wasn't sure what the second item meant in the list of lists.)

I've looked through the github issues associated with the explicit indices, but can't quite tell if I can use them to load FVCOM model output. In any case I just updated and tried without doing anything new and it didn't work:

import xarray as xr url = 'https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/SFBOFS/MODELS/2022/06/30/nos.sfbofs.fields.f014.20220630.t09z.nc' # this file will not be available in a few days but one for the present day will be available ds = xr.open_dataset(url, drop_variables='Itime2')

Same error message as before:

MissingDimensionsError: 'siglay' has more than 1-dimension and the same name as one of its dimensions ('siglay', 'node'). xarray disallows such variables because they conflict with the coordinates used to label dimensions.

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.

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

In order to maintain a list of currently relevant issues, we mark issues as stale after a period of inactivity

If this issue remains relevant, please comment here or remove the stale label; otherwise it will be marked as closed automatically

@rabernat , this unstructured grid model output follows the UGRID Conventions, which layer on top of the CF Conventions. The issue Xarray is having here is with the vertical coordinate however, so this issue could arise with any CF convention model where the vertical stretching function varies over the domain.

As requested, here is the ncdump of this URL: ``` jovyan@jupyter-rsignell-2dusgs:~$ ncdump -h http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc netcdf NECOFS_GOM3_FORECAST { dimensions: time = UNLIMITED ; // (145 currently) maxStrlen64 = 64 ; nele = 99137 ; node = 53087 ; siglay = 40 ; three = 3 ; variables: float lon(node) ; lon:long_name = "nodal longitude" ; lon:standard_name = "longitude" ; lon:units = "degrees_east" ; float lat(node) ; lat:long_name = "nodal latitude" ; lat:standard_name = "latitude" ; lat:units = "degrees_north" ; float xc(nele) ; xc:long_name = "zonal x-coordinate" ; xc:units = "meters" ; float yc(nele) ; yc:long_name = "zonal y-coordinate" ; yc:units = "meters" ; float lonc(nele) ; lonc:long_name = "zonal longitude" ; lonc:standard_name = "longitude" ; lonc:units = "degrees_east" ; float latc(nele) ; latc:long_name = "zonal latitude" ; latc:standard_name = "latitude" ; latc:units = "degrees_north" ; float siglay(siglay, node) ; siglay:long_name = "Sigma Layers" ; siglay:standard_name = "ocean_sigma_coordinate" ; siglay:positive = "up" ; siglay:valid_min = -1. ; siglay:valid_max = 0. ; siglay:formula_terms = "sigma: siglay eta: zeta depth: h" ; float h(node) ; h:long_name = "Bathymetry" ; h:standard_name = "sea_floor_depth_below_geoid" ; h:units = "m" ; h:coordinates = "lat lon" ; h:type = "data" ; h:mesh = "fvcom_mesh" ; h:location = "node" ; int nv(three, nele) ; nv:long_name = "nodes surrounding element" ; nv:cf_role = "face_node_connnectivity" ; nv:start_index = 1 ; float time(time) ; time:long_name = "time" ; time:units = "days since 1858-11-17 00:00:00" ; time:format = "modified julian day (MJD)" ; time:time_zone = "UTC" ; time:standard_name = "time" ; float zeta(time, node) ; zeta:long_name = "Water Surface Elevation" ; zeta:units = "meters" ; zeta:standard_name = "sea_surface_height_above_geoid" ; zeta:coordinates = "time lat lon" ; zeta:type = "data" ; zeta:missing_value = -999. ; zeta:field = "elev, scalar" ; zeta:coverage_content_type = "modelResult" ; zeta:mesh = "fvcom_mesh" ; zeta:location = "node" ; int nbe(three, nele) ; nbe:long_name = "elements surrounding each element" ; float u(time, siglay, nele) ; u:long_name = "Eastward Water Velocity" ; u:units = "meters s-1" ; u:type = "data" ; u:missing_value = -999. ; u:field = "ua, scalar" ; u:coverage_content_type = "modelResult" ; u:standard_name = "eastward_sea_water_velocity" ; u:coordinates = "time siglay latc lonc" ; u:mesh = "fvcom_mesh" ; u:location = "face" ; float v(time, siglay, nele) ; v:long_name = "Northward Water Velocity" ; v:units = "meters s-1" ; v:type = "data" ; v:missing_value = -999. ; v:field = "va, scalar" ; v:coverage_content_type = "modelResult" ; v:standard_name = "northward_sea_water_velocity" ; v:coordinates = "time siglay latc lonc" ; v:mesh = "fvcom_mesh" ; v:location = "face" ; float ww(time, siglay, nele) ; ww:long_name = "Upward Water Velocity" ; ww:units = "meters s-1" ; ww:type = "data" ; ww:coverage_content_type = "modelResult" ; ww:standard_name = "upward_sea_water_velocity" ; ww:coordinates = "time siglay latc lonc" ; ww:mesh = "fvcom_mesh" ; ww:location = "face" ; float ua(time, nele) ; ua:long_name = "Vertically Averaged x-velocity" ; ua:units = "meters s-1" ; ua:type = "data" ; ua:missing_value = -999. ; ua:field = "ua, scalar" ; ua:coverage_content_type = "modelResult" ; ua:standard_name = "barotropic_eastward_sea_water_velocity" ; ua:coordinates = "time latc lonc" ; ua:mesh = "fvcom_mesh" ; ua:location = "face" ; float va(time, nele) ; va:long_name = "Vertically Averaged y-velocity" ; va:units = "meters s-1" ; va:type = "data" ; va:missing_value = -999. ; va:field = "va, scalar" ; va:coverage_content_type = "modelResult" ; va:standard_name = "barotropic_northward_sea_water_velocity" ; va:coordinates = "time latc lonc" ; va:mesh = "fvcom_mesh" ; va:location = "face" ; float temp(time, siglay, node) ; temp:long_name = "temperature" ; temp:standard_name = "sea_water_potential_temperature" ; temp:units = "degrees_C" ; temp:coordinates = "time siglay lat lon" ; temp:type = "data" ; temp:coverage_content_type = "modelResult" ; temp:mesh = "fvcom_mesh" ; temp:location = "node" ; float salinity(time, siglay, node) ; salinity:long_name = "salinity" ; salinity:standard_name = "sea_water_salinity" ; salinity:units = "0.001" ; salinity:coordinates = "time siglay lat lon" ; salinity:type = "data" ; salinity:coverage_content_type = "modelResult" ; salinity:mesh = "fvcom_mesh" ; salinity:location = "node" ; int fvcom_mesh ; fvcom_mesh:cf_role = "mesh_topology" ; fvcom_mesh:topology_dimension = 2 ; fvcom_mesh:node_coordinates = "lon lat" ; fvcom_mesh:face_coordinates = "lonc latc" ; fvcom_mesh:face_node_connectivity = "nv" ;

// global attributes: :title = "NECOFS GOM3 (FVCOM) - Northeast US - Latest Forecast" ; :institution = "School for Marine Science and Technology" ; :source = "FVCOM_3.0" ; :Conventions = "CF-1.0, UGRID-1.0" ; :summary = "Latest forecast from the FVCOM Northeast Coastal Ocean Forecast System using an newer, higher-resolution GOM3 mesh (GOM2 was the preceding mesh)" ;

```