Coordinate computation does not take into accound inverted axis in `center` mode

[Kirill888]

To be compatible with rioxarray one needs to use stackstac.stack(..., xy_coords="center") when computing X/Y coordinate values. When using this mode on data that contains "inverted Y axis", a most common scenario, Y axis coordinates are offset by 1 pixel size into positive direction.

I have made a small reproducer. Data is a global synthetic image with 1 degree per pixel in EPSG:4326, when loading it with xy_coords="center" you would expect Y coordinate to span from -89.5 to 89.5, but instead it goes from 90.5 to -88.5.

https://nbviewer.org/gist/Kirill888/b3dad8afdc10b37cd21af4aea8f417e3/stackstac-xy_coords-error-report.ipynb
https://gist.github.com/Kirill888/b3dad8afdc10b37cd21af4aea8f417e3

This causes issue reported earlier here: #68

Code that performs computation of the coordinate just offsets "top-left" coordinate by positive half pixel size, but instead should offset by sign(coord[1] - coord[0])*abs(resolution)*0.5

[gjoseph92]

@Kirill888 thanks a ton for looking into this. I'm sorry that I haven't had time to work on #68 at all, so I appreciate you digging in.

I wonder if this fixes it:

diff --git a/stackstac/prepare.py b/stackstac/prepare.py
index 8f28ceb..545fcf7 100644
--- a/stackstac/prepare.py
+++ b/stackstac/prepare.py
@@ -400,9 +400,9 @@ def to_coords(
                 half_xpixel, half_ypixel = xres / 2, yres / 2
                 minx, miny, maxx, maxy = (
                     minx + half_xpixel,
-                    miny + half_ypixel,
+                    miny - half_ypixel,
                     maxx + half_xpixel,
-                    maxy + half_ypixel,
+                    maxy - half_ypixel,
                 )
 
             height, width = spec.shape

Reading that code, looks like I forgot that resolutions_xy are always positive; see

stackstac/stackstac/raster_spec.py

Lines 23 to 31 in 3a58906

	def transform(self) -> affine.Affine:
	return affine.Affine(
	self.resolutions_xy[0], # xscale
	0.0,
	self.bounds[0], # xoff
	0.0,
	-self.resolutions_xy[1], # yscale
	self.bounds[3], # yoff
	)

[Kirill888]

@gjoseph92 no worries.

Regarding your change, this will "fix it" for a most common case where Y axis has "negative resolution", i.e. top row has lowest Y coordinate, but it will break for all other cases. GeoTIFF can encode equivalent axis aligned data in 4 different ways depending on which corner is picked as reference. Say you have

xx = xr.DataArray(...)
a = xx[:,:]
b = xx[:,::-1]
c = xx[::-1, :]
d = xx[::-1, ::-1]

Images a,b,c,d will have pixels arranged in four different ways, so if you access them by index they won't match up across all those images, but if you access them by coordinate (assuming center coordinate is used, as this is what xarray expects) you'll see that pixels at the same location have the same value.

So to address this issue for all possible geotiffs out there you really need to compute sign per axis, including X axis, although I have not seen files like that in the wild.

- half_xpixel, half_ypixel = xres / 2, yres / 2
+ sx, sy = numpy.sign(transform.a), numpy.sign(transform.e)
+ half_xpixel, half_ypixel = (s*res/2 for s, res in zip([sx, sy], [xres, yres])

[gjoseph92]

Yeah, I agree that having saying "X resolution is always positive, Y is always negative" is a bit fast and loose. I remember going back and forth on this a bit when I wrote the code initially. I did this because a) resolution, as a positive number, is more intuitive to users to specify and b) when looking at multiple assets at different resolutions and trying to pick the finest, you want them to all be positive numbers (even if some assets are "flipped" and some aren't). Technically, maybe this is more GSD?

I think, though, that in stackstac's case, by construction we really are guaranteed that X is positive and Y is negative. stackstac isn't reading any GeoTIFFs—it's deciding what transform to tell rasterio to use when it reads the GeoTIFFs later. In the implementation of the transform I linked above, we simply decide that e is -self.resolutions_xy[1]—Y is always negative. (And the resolutions_xy, when calculated from STAC metadata, will always be positive.) Eventually, we tell rasterio to warp the input GeoTIFF into this transform and CRS via a VRT, so assuming GDAL handles this correctly, everything should be fine.

So stackstac is always going to warp any orientation of GeoTIFF into positive-X, negative-Y. This may not be what you're expecting, but with that diff I posted, I think it would no longer be wrong: the xarray coordinates would match up to the correct pixels.

Eventually though, we should support respecting the orientation of a transform given by the proj:transform field, so that "flipped" datasets don't get "un-flipped".

[Kirill888]

Thanks for the explanation. I have assumed that .stack returns pixels in native projection/resolution when possible. And native allows for

1. Inverted axis for either X/Y
2. Non-orthographic transforms

So if the desired behaviour is to normalize coordinates to be

1. Rectilinear
2. With inverted Y axis

Then that's fine. But then what are the situations where this code path happens:

stackstac/stackstac/prepare.py

Lines 413 to 420 in 3a58906

	height, width = spec.shape
	if pixel_center:
	xs, _ = transform * (np.arange(width) + 0.5, np.zeros(width) + 0.5)
	_, ys = transform * (np.zeros(height) + 0.5, np.arange(height) + 0.5)
	else:
	xs, _ = transform * (np.arange(width), np.zeros(width))
	_, ys = transform * (np.zeros(height), np.arange(height))

It looks like by construction transform is always rectinlinear with inverted Y axis:

|sx|    0  tx
  0  -|sy| ty
  0     0   1 

stackstac/stackstac/raster_spec.py

Lines 23 to 31 in 3a58906

	def transform(self) -> affine.Affine:
	return affine.Affine(
	self.resolutions_xy[0], # xscale
	0.0,
	self.bounds[0], # xoff
	0.0,
	-self.resolutions_xy[1], # yscale
	self.bounds[3], # yoff
	)

---

I remember going back and forth on this a bit when I wrote the code initially. I did this because a) resolution, as a positive number, is more intuitive to users to specify and b) when looking at multiple assets at different resolutions and trying to pick the finest, you want them to all be positive numbers (even if some assets are "flipped" and some aren't). Technically, maybe this is more GSD?

That one is hard, GDAL/rasterio certainly use negative sign to indicate "direction" when reporting resolution, but I agree that this is awkward in a lot of cases. But I guess this is separate concern from coordinate value computation.

[gjoseph92]

Yeah, I think that entire codepath for if not spec.transform.is_rectilinear could be removed. Because as you point out, the transform that we generate is always rectilinear.

So if the desired behaviour is to normalize coordinates

I don't think it's really the desired behavior... just what the behavior currently is. I plan on keeping that behavior for now and just fixing that issue, but we should better support native projection/resolution in the future.

[gjoseph92]

Closed in #94

[maawoo]

Thanks @gjoseph92 and @Kirill888 for fixing this and #68 !
Also very interesting to follow your discussion.