diff --git a/CHANGELOG.md b/CHANGELOG.md index ae53c3d..ea4b035 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/stackstac/prepare.py b/stackstac/prepare.py index 8f28ceb..b913372 100644 --- a/stackstac/prepare.py +++ b/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 diff --git a/stackstac/raster_spec.py b/stackstac/raster_spec.py index 866266a..2a0f190 100644 --- a/stackstac/raster_spec.py +++ b/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(