Fix one-pixel shift with xy_coords="center"

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## 0.2.2 (...)
 - Support [pystac](https://github.com/stac-utils/pystac) ItemCollections
 - Fix bug where repeated metadata values would be None
+- Fix one-pixel shift when using `xy_coords="center"` [@gjoseph92](https://github.com/gjoseph92) [@Kirill888](https://github.com/Kirill888) [@maawoo](https://github.com/maawoo)
 
 ## 0.2.1 (2021-05-07)
 Support [xarray 0.18](http://xarray.pydata.org/en/stable/whats-new.html#v0-18-0-6-may-2021) and beyond, as well as looser version requirements for some other dependencies.

stackstac/prepare.py

@@ -391,32 +391,27 @@ def to_coords(
             )
 
         transform = spec.transform
-        if transform.is_rectilinear:
-            # faster-path for rectilinear transforms: just `arange` it instead of doing the affine math
-            minx, miny, maxx, maxy = spec.bounds
-            xres, yres = spec.resolutions_xy
-
-            if pixel_center:
-                half_xpixel, half_ypixel = xres / 2, yres / 2
-                minx, miny, maxx, maxy = (
-                    minx + half_xpixel,
-                    miny + half_ypixel,
-                    maxx + half_xpixel,
-                    maxy + half_ypixel,
-                )
+        # We generate the transform ourselves in `RasterSpec`, and it's always constructed to be rectilinear.
+        # Someday, this should not always be the case, in order to support non-rectilinear data without warping.
+        assert (
+            transform.is_rectilinear
+        ), f"Non-rectilinear transform generated: {transform}"
+        minx, miny, maxx, maxy = spec.bounds
+        xres, yres = spec.resolutions_xy
+
+        if pixel_center:
+            half_xpixel, half_ypixel = xres / 2, yres / 2
+            minx, miny, maxx, maxy = (
+                minx + half_xpixel,
+                miny - half_ypixel,
+                maxx + half_xpixel,
+                maxy - half_ypixel,
+            )
 
-            height, width = spec.shape
-            # Wish pandas had an RangeIndex that supported floats...
-            xs = pd.Float64Index(np.linspace(minx, maxx, width, endpoint=False))
-            ys = pd.Float64Index(np.linspace(maxy, miny, height, endpoint=False))
-        else:
-            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))
+        height, width = spec.shape
+        # Wish pandas had an RangeIndex that supported floats...
+        xs = pd.Float64Index(np.linspace(minx, maxx, width, endpoint=False))
+        ys = pd.Float64Index(np.linspace(maxy, miny, height, endpoint=False))
 
         coords["x"] = xs
         coords["y"] = ys

stackstac/raster_spec.py

@@ -19,6 +19,15 @@ class RasterSpec:
     bounds: Bbox
     resolutions_xy: Resolutions
 
+    def __post_init__(self):
+        xres, yres = self.resolutions_xy
+        assert xres > 0, f"X resolution {xres} must be > 0"
+        assert yres > 0, f"Y resolution {yres} must be > 0"
+
+        minx, miny, maxx, maxy = self.bounds
+        assert minx < maxx, f"Invalid bounds: {minx=} >= {maxx=}"
+        assert miny < maxy, f"Invalid bounds: {miny=} >= {maxy=}"
+
     @cached_property
     def transform(self) -> affine.Affine:
         return affine.Affine(