Adressing floating-point issues

[benbovy]

I'm reading again the pangeo thread about the limitations of NetCDF-style coordinates for large geospatial rasters to see whether RasterIndex in its current implementation would already fix some of the problems outlined in that thread, including floating-point issues.

This comment makes me think that "just" keeping track of the affine parameters within RasterIndex will probably not solve all those floating-point issues.

For example, when extracting a Region Of Interest (ROI) by data slicing, RasterIndex currently computes new affine parameters that are then wrapped in a new RasterIndex instance. However, if the affine parameters are float numbers, this will not ensure that the pixels of the extracted ROI will be perfectly aligned with those of the original dataset (I think?).

In the comment linked above, there are two suggestions for addressing this problem:

1. Keep track of both the original affine parameters and slice offsets instead of computing new affine parameters
2. Introduce some alignment tolerance

I think either suggestion could be pretty easily supported by RasterIndex (I don't see any big technical obstacle for implementing them).

[dcherian]

Yes, I agree. Though I think the problem is much better in that we now only compare the "offset" (or origin) so two points instead of 1000s of points.

[dcherian]

Here is a concrete example. This dataset, GLAD Land Use Land Cover is distributed as tiles.

So a user should be able to read these 4 TIFFs and assemble them in a 2x2 grid to form a combined dataset.

https://storage.googleapis.com/earthenginepartners-hansen/GLCLU2000-2020/v2/2000/00N_040W.tif
https://storage.googleapis.com/earthenginepartners-hansen/GLCLU2000-2020/v2/2000/00N_050W.tif
https://storage.googleapis.com/earthenginepartners-hansen/GLCLU2000-2020/v2/2000/10S_050W.tif
https://storage.googleapis.com/earthenginepartners-hansen/GLCLU2000-2020/v2/2000/10S_050W.tif

---

Example

import rioxarray
import xarray as xr

rioxarray.set_options(use_raster_index=True)
ds40 = xr.open_dataset("glad/00N_040W.tif", engine="rasterio")
ds50 = xr.open_dataset("glad/00N_050W.tif", engine="rasterio")
xr.concat([ds40, ds50], dim=("x", "y"))

Interestingly I have to specify dim=("x", "y") to make this work. Simply dim="x" raises a hard-to-interpret error:

ValueError: cannot exclude dimension(s) x from alignment because these are used by an index together with non-excluded dimensions y

• should we opt the user in to ("x", "y") if either "x" or "y" are provided?
• or simply have the error message suggest that the user pass in ("x", "y")

Note that this is the really simple case. The two Affines are:

Affine(0.00025, 0.0, -40.0,
       0.0, -0.00025, 0.0)

Affine(0.00025, 0.0, -50.0,
       0.0, -0.00025, 0.0)