How to recognise a rectangular grid?

[ChrisBarker-NOAA]

The current code for recognising a rectangular grid is simple, but broken:

    def recognize(ds: xr.Dataset) -> bool:
        """Recognize if the dataset matches the given grid."""
        lat = ds.cf.coordinates.get("latitude", None)
        lon = ds.cf.coordinates.get("longitude", None)
        if lat is None or lon is None:
            return False

        # Make sure the coordinates are 1D and match
        lat_ndim = ds[lat[0]].ndim
        lon_ndim = ds[lon[0]].ndim
        return lat_ndim == lon_ndim and lon_ndim == 1

All this does is:

• make sure there are a latitude and longitude coordinate -- fair enough
• make sure that the first of such are 1-D -- god start, but ....

If you run this code in an unstructured grid, it returns True -- not good.

Why?

Because UGRids have a bunch of 1-D lat and long coords:

In [18]: ds.cf.coordinates
Out[18]: 
{'longitude': ['mesh_boundary_lon',
  'mesh_edge_lon',
  'mesh_face_lon',
  'mesh_node_lon'],
 'latitude': ['mesh_boundary_lat',
  'mesh_edge_lat',
  'mesh_face_lat',
  'mesh_node_lat']}

and they are all 1D.

So: what do I test for?

1. this code looks at the zeroth lat-lon only -- is that OK? in my examples there are only one of each, but would that always be the case?

2. if above, it seem we should look and see if variables are using both the lat and lon coordinates ...

But this does seem fragile ...

Anyone have a better idea?

@ocefpaf: any thoughts on this?

@davidhassell: Sorry to ling you on this, but I thught you might have some ideas.

[ChrisBarker-NOAA]

OK -- worked out a solution that seems to work:

    def recognize(ds: xr.Dataset) -> bool:
        """
        Recognize if the dataset matches the given grid.
        """
        lat = ds.cf.coordinates.get("latitude", None)
        lon = ds.cf.coordinates.get("longitude", None)
        if lat is None or lon is None:
            return False

        # choose first one -- valid assumption??
        lat = lat[0]
        lon = lon[0]
        # Make sure the coordinates are 1D and match
        if not (1 == ds[lat].ndim == ds[lon].ndim):
            return False

        # make sure that at least one variable is using both the
        #   latitude and longitude dimensions
        #   (ugrids have both coordinates, but not both dimensions)
        for var_name, var in ds.data_vars.items():
            if (lon in var.dims) and (lat in var.dims):
                return True
        return False

What makes me most nervous is arbitrarily choosing the first latitude and longitude coordinate found. So far, the rectangular grid examples I've found only have one, but ???

[ChrisBarker-NOAA]

This is in branch: regular_grid_tests. I'm hoping to merge that in soon, but any ideas better are welcome.

I also added more tests, so that the three grid types all make sure that the recognise themselves, and don't recognize the others ... at least with the small number of test files we have.

[Nesar976]

Hi @ChrisBarker-NOAA
I feel the check should be based more on how the data is shaped rather than just the presence of 1-D lat/lon.

If there is at least one variable that actually uses both latitude and longitude as dimensions (for example var(lat, lon)), then it makes sense to treat it as a rectangular grid. If lat/lon only describe points or mesh elements and are never used together, it probably shouldn’t be classified as rectangular.

This would avoid arbitrarily picking the first lat/lon and should also prevent unstructured grids from being misidentified.

[ocefpaf]

@ChrisBarker-NOAA I don't believe we have enough metadata to do this unambiguously, right? Or at least as a first check and then, if absent, fallback to that one above? I have no experience with triangular grids, but I'll try to gather some netcdf files (I believe we have some in the ugrid repos?) and experiment a bit.

[ocefpaf]

Chris, one crazy idea would be to try to create the triang object from matplotlib, or something similar, with a subset of the grid to avoid making this check super expensive.

[ChrisBarker-NOAA]

If there is at least one variable that actually uses both latitude and longitude as dimensions (for example var(lat, lon)), then it makes sense to treat it as a rectangular grid.

That's what my current version does -- see the code above.

This would avoid arbitrarily picking the first lat/lon

well, to see if there are anyvariables that use both, you need to know what dimension to use, so you need to have a lat and lon to use -- my version so far uses the first one of each (and all the test files only have one, so that's easy).

I suppose we could do a loop and check all the combinations of all the latitude and longitude variables -- a bit ugly, but it would be robust. PRs welcome if someone wants to write that :-)

(though make sure you have a test case that is rectangular and has more than one latitude and longitude coordinate)

Current draft in the rectangular_grid_tests branch -- which I may merge soon -- at least it's better than what's there now, which is quite broken.

https://github.com/ioos/xarray-subset-grid/blob/regular_grid_tests/xarray_subset_grid/grids/regular_grid.py#L96

Chris, one crazy idea would be to try to create the triang object from matplotlib, or something similar, with a subset of the grid to avoid making this check super expensive.

well, that would be a way to check if it's a traingular mesh, but not if it's something else that isn't rectangular.

If we did something liike that, it would be in the ugrid recogniser. -- that one currently looks for the UGRID compiliant mesh variable. It might be nice to have some auto-detection, but now a user can call ugrid.assign_ugrid_topology first and it will work fine.

[davidhassell]

I'm possibly missing a trick, but can you check for the existence of the mesh attribute on the data variable (or perhaps xarray doesn't save it?)

[ocefpaf]

or perhaps xarray doesn't save it?

It should be under ds.attrs, but the problem is how many models out there define this?

[ChrisBarker-NOAA]

I'm possibly missing a trick, but can you check for the existence of the mesh attribute on the data variable

that's exactly what the UGRID detector does.

But in this case, this is the rectangular grid detector -- it shouldn't return a True if it's not a valid rectangular grid, even if it's not a valid UGRID .... it could be who knows what else?

In fact, that's the use-case that I found this bug with -- I had a UGRID file that didn't have the mesh defined -- and it got incorrectly recognised as a rectangular grid -- which then barfed when you tried to use it.

For a bit of background, this project provides a way to subset various gridded model results -- it need to know what sort of grid it is so that it can be subset properly. So there are different objects for each grid type, and each one has a "recognise" method. Ideally, it will return True is the dataset is its grid type, and false if not. The code loops through all the grid types until it finds one that recognises it.

@davidhassell -- I pinged you because I thought you might have an idea of the clear definition for a rectangular grid in CF.

Would that be if there was at least one variable with a (longitude, latitude) in it's dimensions?

The question is: there can be multiple latitude and longitude coordinates -- so should we chack againsta ll of them for a variable with one of each?

I have code that works for all the examples we have, but that's not a comprehensive list at all -- and Ideally I'd like to get something robust in while I'm thinking about it.

[davidhassell]

Oh , I see, now! I'll have a think (def recognize looks pretty good, though)