Add some notes on RasterIndex design

[dcherian]

cc @benbovy @maxrjones @scottyhq

[benbovy]

Another complication is the special but very common case where the model space can be represented by orthogonal 1D coordinates (e.g., rectilinear affine with no rotation), which by design is not supported by the xarray.CoordinateTransform base class.

In corteva/rioxarray#846 the workaround consists in introducing an AxisAffineTransform(CoordinateTransform) subclass that wraps a 2D affine but handles one axis, as well as a AxisAffineTransformIndex(CoordinateTransformIndex) subclass. The latter allows providing useful features in the 1D case, such as automatically converting the index into a PandasIndex when indexing at arbitrary locations.

[benbovy]

I'm not familiar with RPCs, but I wonder if this wouldn't allow an Xarray Dataset with, e.g.,

• x and y coordinates associated with a RasterIndex[AffineTransformIndex]

and

• lat and lon coordinates associated with a RasterIndex[RPCIndex].

which would enable both ds.sel(x=..., y=...) and ds.sel(lat=..., lon=...).

However, I read that

The RPC model allows a row/column location to be computed for a given lat, long and height value. It is not inherently invertable, though it is usually possible to compute at lat,long location from a row, column and height value using iterative methods.

So IIUC CoordinateTransform.forward() would be trickier to implement. Maybe we could make the implementation optional? This would mean that lat, lon are virtual coordinates that cannot be materialized. But since the RPC model allows CoordinateTransform.inverse() it would still allow ds.sel(lat=..., lon=...).

[dcherian]

This would mean that lat, lon are virtual coordinates that cannot be materialized.

Yes, that's right. I think this may still work with the current model we have, perhaps we need to override the repr.

[benbovy]

Is this class meant to inherit from xarray.CoordinateTransform or is it a different concept?

[dcherian]

It's a different concept that represents which of the 5 ways the transformation is recorded in GeoTIFF metadata.

[benbovy]

That looks great @dcherian. I left a few comments mainly about how to build RasterIndex on top of Xarray coordinate transform, although I'm not sure that these are really relevant here since you seem to focus the notes on how to pass information from/to metadata.

[benbovy]

Apparently a 6th possible one could be based on subsampled auxiliary coordinates, detailed in CF section 8.3 and equivalent to GDAL's geolocation arrays with PIXEL_STEP and/or LINE_STEP > 1.

[benbovy]

I think that a decode/encode workflow strategy makes a lot of sense in this case.

Taking the example 8.3 from CF section 8.3, the decode step may consist in:

• turn the tie point coordinate variables lat(tp_yc, tp_xc) and lon(tp_yc, tp_xc) into lat(yc, xc) and lon(yc, xc) Xarray coordinates associated with a custom transformation index that stores only the tie points. In other words, uncompress the dimensions of the lat & lon coordinates without uncompressing their data.
• also remove the tie point index variables and the interpolation variable, and track their data / metadata internally in the index

The encode step would then consist in restoring the compressed tie point coordinate & index variables as well as the interpolation variable.

[dcherian]

Ah nice, i agree completely