Implement cell bounds convention for xarray

In GitLab by @Huite on May 18, 2021, 14:34

It looks like xarray finally has some support for cell boundaries arrays: https://github.com/pydata/xarray/pull/2844

The default option is simply storing cell "midpoints", but this is strictly insufficient to derive the boundaries, if cell size is not a constant (like it is in general GIS rasters).

MODFLOW supports non-equidistant grids, and thus IDFs (and of course netCDFs) can be non-equidistant. To support regridding in these cases, I had to invent my own convention of storing d coordinates for the cell sizes (dz, dy, dx). Unfortunately, while there is some precedent for scalar cell sizes (e.g. from xr.open_rasterio), this is not the case for array cell sizes.

So rather than storing dy and dx arrays for non-equidistant data, we could simply store the cell boundaries instead. Note that cell centers are still required as a coordinate, since software like GDAL or MDAL will not look at boundaries at all. In some cases, dx and dy are useful (e.g. to compute cell area), but having an ad hoc / informal convention like this is a big downside.

I guess we could slowly deprecate dx and dy (and dz) coordinate arrays?

This has no relation with either MDAL or GDAL: currently they do not understand either. Given the fact, however, that cell bounds is CF compliant, while dx and dy coordinates are not, makes it a lot more likely for bounds support to be included at some point in the future.

Some changes:

• Convert IDF dx and dy arrays to cell bounds
• Maintain scalar dx and dy, since xr.open_rasterio still does this (also means that most use will not be impacted, since non-equidistant grids are quite rare)
• Deprecate (in due time) the dx and dy array based logic. This removes the duck typing complexities for dx and dy, since they can always be assumed to be scalar.
• For the regridder, it no longer needs to infer cell boundaries, they are assumed given. For scalar dx, they are (still) easily computed. Will still need an equidistance check (but checks for presence of cell boundaries, rather than presence of dx dy coords).
• For IDF output: automatically convert cell boundary coordinates to the required dx and dy arrays in case of non-equidistance.

Start by tracing the IDF import and output paths and the regridder's handling of xarray coordinates. Review the CF cell-boundary convention and the existing dx, dy, and dz behavior. Done means non-equidistant grids use cell boundaries consistently, scalar spacing remains supported, and the listed IDF and regridding cases work without relying on array-based spacing.