Example which highlights the limitations of NetCDF-style coordinates for large geospatial rasters

This issue comes up every now and then:

It’s important to understand that GeoBox is recomputed from coordinates, that’s needed to support slicing into geo-referenced data, and also to support data constructed by other libraries. BUT there is no guarantee that this recomputed GeoBox will produce exactly the same coordinates when used to create a new array from it. Essentially GeoBox -> coords -> GeoBox is not guaranteed to be lossless. It can only be lossless when both resolution and translation components of the Affine matrix are basically integers. Sentinel-2 has that property for example, scale=+/-10 tx,ty=10*N where N is an integer.

Not sure what a proper solution for this should be. We can keep track of the original resolution in an attribute of the coordinate, and use exactly that value when extracting GeoBox from coords (with a check to deal with xx[::10, ::10] type of slicing. I guess we can also keep track of original translation, and only recompute that if array has been sliced. I’ll probably implement that for the next version of odc-geo, actually.

The problem of creating sub-geobox that will be able to produce exactly the same coords as the original geobox in the sliced section is much harder to address. We want this invariant:

gbox[roi].coords == gbox.coords[roi]

That’s not possible unless gbox[roi] retains parent and slice, and essentially return self.parent.coords[self.roi]. Or we implement “rounding to some fraction of a pixel” kinda logic.