Can't subset RTOFS - doesn't recognize grid

[AmyMacFadyen]

What went wrong?

It throws this warning
UserWarning: no grid type found in this dataset
warnings.warn("no grid type found in this dataset")

And returns None

Also of note -- you only get the warning if you open a multifile dataset -- but it returns None regardless (i.e., it doesn't work in either case)

Operating System

Windows

Version

0.1.1.dev8+g97ec48035

Python Version

3.10

Code to Reproduce

ds = xr.open_mfdataset(['http://edsdata.oceansmap.com/thredds/dodsC/EDS/HYCOM_GLOBAL/hycomglobal_20260619.nc','http://edsdata.oceansmap.com/thredds/dodsC/EDS/HYCOM_GLOBAL/hycomglobal_20260620.nc'])
bbox = (-126.0, 47.5, -124.0, 48.5)
ds.xsg.subset_bbox(bbox)

Errors, Traceback, and Logs

In [106]: ds.xsg.subset_bbox(bbox)
C:\Users\amy.macfadyen\PyProjects\xarray-subset-grid\xarray_subset_grid\accessor.py:47: UserWarning: no grid type found in this dataset
  warnings.warn("no grid type found in this dataset")

In [107]:

[ChrisBarker-NOAA]

To make it a bit easier, here's the header:

dimensions:
	time = UNLIMITED ; // (9 currently)
	depth = 1 ;
	x = 4500 ;
	y = 3298 ;
variables:
	float time(time) ;
		time:units = "days since 1900-12-31 00:00:00" ;
		time:calendar = "standard" ;
		time:_ChunkSizes = 1024 ;
	short water_u(time, depth, y, x) ;
		water_u:long_name = "Eastward Water Velocity" ;
		water_u:add_offset = "0.f" ;
		water_u:scale_factor = 0.001f ;
		water_u:_FillValue = -30000s ;
		water_u:units = "m/s" ;
		water_u:_ChunkSizes = 1, 1, 1649, 2250 ;
	short water_v(time, depth, y, x) ;
		water_v:long_name = "Northward Water Velocity" ;
		water_v:add_offset = "0.f" ;
		water_v:scale_factor = 0.001f ;
		water_v:_FillValue = -30000s ;
		water_v:units = "m/s" ;
		water_v:_ChunkSizes = 1, 1, 1649, 2250 ;
	float lat(y, x) ;
		lat:long_name = "Latitude" ;
		lat:units = "degrees_north" ;
		lat:axis = "Y" ;
	float lon(y, x) ;
		lon:long_name = "Longitude" ;
		lon:units = "degrees_east" ;
		lon:axis = "X" ;
	float depth(depth) ;
		depth:long_name = "Depth" ;
		depth:units = "meter" ;
		depth:axis = "Z" ;

So we have a 2D specification for a logically rectangular grid -- which I think should be readable as a "SGRID".

I suspect it doesn't detect it because the dimension names are 'x' and 'y' rather than 'lat' and 'lon' or some such.

But we'll need to see what the grid detection actually does here.

I think the first step is to make it CF compliant, which I don't think it is. And then the current code may work.

@ocefpaf: Is this something you could take a look at?

[ocefpaf]

Sure. What would be the goal of xarray-subset-grid when the issue is a non-compliant file?

1. Document on how to patch the at loading time. (I need to check if that is possible here, I do know that both xarray and iris can do some metadata patching);
2. Ask users to fix it with a CF/compliance error message;
3. Handle it internally?

I'm not a big fan of the last one. We could end up with tons of code just to work around data. I'll see if I can do 1 easily later today and post an example.

---

Edit: The grid seems to be 3D, not 2D

<xarray.DataArray 'lon' (time: 18, y: 3298, x: 4500)> Size: 1GB
dask.array<concatenate, shape=(18, 3298, 4500), dtype=float32, chunksize=(9, 3298, 4500), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 144B 2026-06-19 ... 2026-06-21
Dimensions without coordinates: y, x
Attributes:
    long_name:  Longitude
    units:      degrees_east
    axis:       X

[ocefpaf]

Here is a quick assessment of what I understand so far (sorry, I'm still catching up with xarray-subset-grid code):

import cf_xarray
import xarray_subset_grid
import xarray as xr


fnames = [
    "http://edsdata.oceansmap.com/thredds/dodsC/EDS/HYCOM_GLOBAL/hycomglobal_20260619.nc",
    "http://edsdata.oceansmap.com/thredds/dodsC/EDS/HYCOM_GLOBAL/hycomglobal_20260620.nc",
]

ds = xr.open_mfdataset(fnames)
ds  # check that the coordinates are only `time` and `depth`.

ds.cf["X"]  # this will fail b/c the data is not CF compliant.

#  can assign coordinates to the x and x dimensions with:

ds = ds.assign_coords(x=ds["lon"])
ds = ds.assign_coords(y=ds["lat"])
ds  # the problem here is that `x and `y` are `(time, y, x)`!
ds.cf["X"]  # but this works now!

However, I still get a no grid found when trying to use subset_bbox:

bbox = (-126.0, 47.5, -124.0, 48.5)
ds.xsg.subset_bbox(bbox)  # This should be an error instead of returning None and a warning BTW

---

Edit: RegularGrid expects 1D and SGrid fails b/c there is no grid_topology (ds.cf.cf_roles["grid_topology"]).

[ChrisBarker-NOAA]

What would be the goal of xarray-subset-grid when the issue is a non-compliant file?

Document on how to patch the at loading time.

This one.

The goal is for xsg to work with compliant datasets. So it's up to the user to ensure that the datasets are compliant first.

It does provide some utilities to make that easier - but users need to call that explicitly.

E.g. ugrid.assign_topology()

[ChrisBarker-NOAA]

RegularGrid expects 1D and SGrid fails b/c there is no grid_topology (ds.cf.cf_roles["grid_topology"]).

First -- this, I think, was created from grib files via TDS. So the first thing we should look at is whether it really needs a 2D lat-lon grid spec. It may not. In which case the fix is to replace the 2D lat-lon with 1D lat lons.

If it really does need 2D grid spec, then it needs a grid_topology variable -- maybe time t0 start a sgrid.assign_topology function.

Also: there is this issue: #104 -- which is not fully resolved.

Thanks!

[AmyMacFadyen]

One thing I realized after posting this -- when I used open_mfdataset, it is concatenating variables that are not time dependent along the time dimension also. That's the default behaviour in xarray which I'm not a fan of. So the lon/lat variables are 3d in that case. I would recommend troubleshooting the case using a single file with xarray.open_dataset

[AmyMacFadyen]

And the grid is curvilinear north of 47N so 2d lat/lon is needed

[ChrisBarker-NOAA]

Thanks Amy -- so I think we need a sgrid.assign_topology() utility function :-)

The MFdataset aggregation is another issue, that I think is outside of XSG's purview, but documenting how to address it would be great.

-CHB

[ocefpaf]

We can make .assign_topology() be something that attaches the metadata. It could be something like:

attrs = {
    "cf_role": "grid_topology",
    "topology_dimension": 2,
    "node_dimensions": "y x",
    "face_dimensions": "",
}

def  assign_topology(self, attrs):
    grid = xr.DataArray(
        data=1,
        name="grid",
        attrs=attrs,
    )
    self.ds["grid"]  = grid
    return

I can't experiment with that at the moment b/c the data is not accessible for me when I try to compute anything on it. This can be check if you try to fully load it and save:

url = "http://edsdata.oceansmap.com/thredds/dodsC/EDS/HYCOM_GLOBAL/hycomglobal_20260619.nc"

ds = xr.open_dataset(url)
ds.to_netcdf("tmp.nc")

That errors out with ERR: DAP DATADDS packet is apparently too short.

@AmyMacFadyen can you access the data using this code block? If so, can you send me a single file for testing?

[ChrisBarker-NOAA]

I'll send you a single file -- not using MFDataset -- nccopy worked for me.

• GitHub doesn't seem to be allowing me to attach it here.

NOTE: it should maybe go in out test files collection, but I've lost track of how we're handling that.

[ChrisBarker-NOAA]

@ocefpaf: take a look at the ugrid.assign_topology -- it builds in the UGRID structure, rather than requiring users to know what to do -- I think that's a good approach.

The first version doesn't have to be full featured though ;-).

[ocefpaf]

@ChrisBarker-NOAA I'm not an expert on SGRID, so apologies if this question doesn't make sense, but when reading the code I found:

xarray-subset-grid/xarray_subset_grid/grids/sgrid.py

Lines 178 to 180 in 709b7b9

	coord_str = grid_topology.attrs.get(f"{location}_coordinates", None)
	if dim_str is None or coord_str is None:
	raise ValueError(f"Could not find {location} dimensions or coordinates")

and that would lead to an error in parsing b/c coord_str is not found. However, that information is optional in the SGRID definition, right?

[ocefpaf]

BTW, I wonder if we should try to guess the grid type when slicing unknown models or if we should default to a basic xarray slice:

ds.where(
    (xmin <= ds["lon"]) & (ds["lon"] <= xmax) &
    (ymin <= ds["lat"]) & (ds["lat"] <= ymax),
    drop=True,
)

The logic would be to implement .sel and .where for 1D and 2D respectively when grid is unknown and issue a warning.

[ChrisBarker-NOAA]

First: I didn't write that code, nor look at it much, so I'm as blind as you are :-) But I think this code was written more by following an example file(s) than the spec ...

However, that information is optional

yes, it's optional, but I think it should not be optional if the coordinates do exist -- that is, having teh coordinates is optional, not having them named if they are there.

And in this case, we need to know the coordinates, so that we can include them in the subset.

But what should the error be?

If the coordinates aren't there at all, then we can't subset them anyway, so maybe this isn't an error, but rather a warning? Maybe something like:

If dim_str is None or coord_str is None: 
    warning.warn(f"Dimension or coordinates for {location} not found. It should be specified in the `{location}_coordinates` attribute of the `grid_topology` variable").

[ocefpaf]

If the coordinates aren't there at all, then we can't subset them anyway, so maybe this isn't an error, but rather a warning? Maybe something like:

Makes sense. I'm trying to wrap my head around what would be the proper metadata patch for this particular case though. I'll get back to this issue ASAP.