Unify 1D and 2D unknown grids and error out when no grid is present

[ocefpaf]

@ChrisBarker-NOAA this PR implements:

• Unified subset for "generic" unknown 1D and 2D grids for datasets that are CF compliant (have at least latitude and longitude recognized by cf-xarray);
• Error out when a grid is not found instead of returning None.

What are we missing?

• We need to rename a few more classes and function to unknown 1D/2D grids instead of regular;
• Add a proper testing for the unknown 2D bbox case.
• Add a test for the polygon case (should be NotImplemented for now.

I did use the rtofs_example.nc file for the test b/c there is something odd with the data variables. However, the logic in this PR works and one can test it with:

import cf_xarray
import xarray_subset_grid
import xarray as xr

ds = xr.open_dataset("rtofs_example.nc")

# This is the bare minimum metadata so cf_xarray can recognize the axis.
ds["lon"] = (((ds["lon"] + 180) % 360) - 180)
ds = ds.assign_coords(lon=ds["lon"])
ds = ds.assign_coords(lat=ds["lat"])
ds["lon"].attrs.update({"standard_name": "longitude"})
ds["lat"].attrs.update({"standard_name": "latitude"})

# Works as expected.
ds[["lon"]].xsg.subset_bbox(bbox) 

# Fails due to a problem with the dataset.
ds[["water_u"]].xsg.subset_bbox(bbox) or ds.xsg.subset_bbox(bbox)

# One can see this when trying to realize the data in that variable:
ds["water_u"].to_numpy()

Partially adress #144
Closes #104

[ocefpaf]

We should probably rename all the classes and functions that refer to "regular grid" to be unknown (1D/2D) grid b/c the logic to subset them is the same. However, that is a breaking change and I want to be sure all devs are onboard before we rename everything.

[ChrisBarker-NOAA]

I don't like "unknown" -- but "regular grid" is a misnomer -- so I'm all for renaming it "rectangular grid".

And it it's a 2D lat-lon, then it's not a rectangular grid, it's curvilinear, which is what SGRID handles, though it's the simplest case of SGRID. It might be hard to have an SGRID-with-a-proper grid spec and simple curvilinear grid be the same class -- maybe a mix-in for the curvilinear grid logic?

[ocefpaf]

It might be hard to have an SGRID-with-a-proper grid spec and simple curvilinear grid be the same class -- maybe a mix-in for the curvilinear grid logic?

I like the "special case inside SGRID". It is SGRID after all and, even if we crack the fix-the-metadata route, the code to subset would be the same b/c it is likely the most accurate there (TBD).

[ChrisBarker-NOAA]

well, the trick is that, I think, the logic references the grid variable to know what's what. So we need to either create a grid variable (which I don't like doing auto-magically, but maybe ??) or there needs to be different logic.

It is SGRID after all

well, kinda -- if there's nothing staggered, then it's nota staggard grid -- it is a specail case of an SGRID, but since the SGRID spec isn't in CF (yet), and it's a perfectly valid CF grid, then we should handle it somewhat natively.

personally, I like to use subclassing to handle "special cases", but there's multiple ways to go there.

[ocefpaf]

We should look at how rioxarray does their polygon subsets. Maybe we can use that instead of custom code.

[ChrisBarker-NOAA]

Worth a look, but oddly the docs don't seem to describe what rioxarray is at a very high level, but it looks like it's for rasters, which is not what we need here :-(

The big distinction is that a raster is a truly regular grid -- it can be defined by a corner, number of pixels and a delta x and delta y.

whereas most of our use cases have delta x and/or delta y that vary.

And, oddly, CF doesn't have a way to support a "proper" raster at all.

And we need the custom polygon code for curvilinear and unstructured grids anyway.

Though now that I think about it -- if you are subsetting a north-aligned rectangular grid, and returning a rectangular grid, then all you can do is make the subset the bounding box of the polygon anyway.

[ocefpaf]

Worth a look, but oddly the docs don't seem to describe what rioxarray is at a very high level, but it looks like it's for rasters, which is not what we need here :-(

Sorry, should've been more specific. It is focused on sub-setting raster images loaded as xarray datasets. The subset method (clip in their method) is generic enough to work with CF datasets that have a project. Ideally this method should be a stand alone library that both rioxarray and xarray-subet-grid could use. However, it is likely that we will re-create their methods centered at CF logic.

[ocefpaf]

@ChrisBarker-NOAA here is an example on how to use a polygon to subset with rioxarray. If it wasn't for the rasterio dependency, b/c t brings in gdal, we could just depend on rioxarray. That is why I prefer to re-write their sub-setting mechanism here.

https://gist.github.com/ocefpaf/e0b56132d3ad1d81f020dd18062d4d21

PS: I'm not sure I used xarray-subset-grid polygon feature correctly there. I cannot get the same triangle as rioxarray's clip.

[ChrisBarker-NOAA]

NOTE: we should finish #104 -- I kind of left that dangling, but I think it's ready to go.

[ocefpaf]

NOTE: we should finish #104 -- I kind of left that dangling, but I think it's ready to go.

Indeed. I missed that, but this implements the same logic to avoid false positives with ugrid.

[ChrisBarker-NOAA]

I think rioxarray's clip is a different idea than subset_grid's subset.

We are not trying to mask out what's not in the polygon -- we're trying to get a full rectangular subset that includes the entire polygon.

In truth subsetting to a polygon is not very useful for rectangular grids -- what you get is the bounding box of the polygon -- which is pretty simple.

The real use-case for polygon subsetting is unstructured grids -- and curvilinear grids in between. Providing it for rectangular grids is really just to get a single API.

[ocefpaf]

The real use-case for polygon subsetting is unstructured grids -- and curvilinear grids in between. Providing it for rectangular grids is really just to get a single API.

In that case I believe that the method with the .where would suffice, right?

[ChrisBarker-NOAA]

In that case I believe that the method with the .where would suffice, right?

Probably, yes -- we just need to make sure to get the edges right.

And something sane for a tiny subset (one grid cell or less ...)

[ocefpaf]

Probably, yes -- we just need to make sure to get the edges right.

I'll make a few experiments to see if we are doing things correctly there.

And something sane for a tiny subset (one grid cell or less ...)

That is challenging. From a pure subset point of view we should return 1 cell if there is any overlap. However, are in a territory where some sort of 2D interpolation would yield more accurate results. With that said, the former seems, IMO, to fall within the goals of this library while the latter does not.

---

Edit: I'm making this one a draft again to wait on #140 and avoid major conflicts later.

[ChrisBarker-NOAA]

With that said, the former seems, IMO, to fall within the goals of this library while the latter does not.

agreed -- the goal here is to get a "Minimal" subset, but to preserve the original grid, interpolation / re-gridding is a whole other topic.

Is a one-cell grid valid? I have to think about that more, if not, then maybe the smallest returned grid could be a 2 x 2 (for structured grids, anyway).