Fix ascending-y grids: affine from coords, and rasterize/mask/coverage engines

src/rasterix/rasterize/core.py

@@ -234,6 +234,20 @@ def rasterize(
         obj = clip_to_bbox(obj, geometries, xdim=xdim, ydim=ydim)
 
     affine = _get_affine(obj, x_dim=xdim, y_dim=ydim)
+    if affine.e > 0:
+        # Engines assume north-up; rasterize the flipped grid and flip back.
+        flip = {ydim: slice(None, None, -1)}
+        return rasterize(
+            obj.isel(flip),
+            geometries,
+            engine=engine,
+            xdim=xdim,
+            ydim=ydim,
+            all_touched=all_touched,
+            merge_alg=merge_alg,
+            geoms_rechunk_size=geoms_rechunk_size,
+            **engine_kwargs,
+        ).isel(flip)
     engine_merge_alg = _normalize_merge_alg(merge_alg, resolved_engine)
 
     rasterize_geometries, dask_rasterize_wrapper = _get_rasterize_funcs(resolved_engine)
@@ -382,6 +396,20 @@ def geometry_mask(
         obj = clip_to_bbox(obj, geometries, xdim=xdim, ydim=ydim)
 
     affine = _get_affine(obj, x_dim=xdim, y_dim=ydim)
+    if affine.e > 0:
+        # Engines assume north-up; mask the flipped grid and flip back.
+        flip = {ydim: slice(None, None, -1)}
+        return geometry_mask(
+            obj.isel(flip),
+            geometries,
+            engine=engine,
+            xdim=xdim,
+            ydim=ydim,
+            all_touched=all_touched,
+            invert=invert,
+            geoms_rechunk_size=geoms_rechunk_size,
+            **engine_kwargs,
+        ).isel(flip)
 
     np_geometry_mask, dask_mask_wrapper = _get_mask_funcs(resolved_engine)
 

src/rasterix/rasterize/exact.py

@@ -308,6 +308,17 @@ def coverage(
         obj = clip_to_bbox(obj, geometries, xdim=xdim, ydim=ydim)
 
     affine = _get_affine(obj, x_dim=xdim, y_dim=ydim)
+    if affine.e > 0:
+        # exactextract assumes north-up; compute on the flipped grid and flip back.
+        flip = {ydim: slice(None, None, -1)}
+        return coverage(
+            obj.isel(flip),
+            geometries,
+            xdim=xdim,
+            ydim=ydim,
+            strategy=strategy,
+            coverage_weight=coverage_weight,
+        ).isel(flip)
     shape = (obj.sizes[ydim], obj.sizes[xdim])
 
     if is_in_memory(obj=obj, geometries=geometries):
@@ -370,8 +381,7 @@ def coverage(
         for c in xy_coords:
             coords = coords.merge(c)
         coords = coords.coords
-    coverage = xr.DataArray(dims=("geometry", ydim, xdim), data=out, coords=coords, attrs=attrs, name=name)
-    return coverage
+    return xr.DataArray(dims=("geometry", ydim, xdim), data=out, coords=coords, attrs=attrs, name=name)
 
 
 # ============================================================================

src/rasterix/utils.py

@@ -147,9 +147,7 @@ def get_affine(
 
     dx = (x[1] - x[0]).item()
     dy = (y[1] - y[0]).item()
-    return Affine.translation(
-        x[0].item() - dx / 2, (y[0] if dy < 0 else y[-1]).item() - dy / 2
-    ) * Affine.scale(dx, dy)
+    return Affine.translation(x[0].item() - dx / 2, y[0].item() - dy / 2) * Affine.scale(dx, dy)
 
 
 def get_crs_from_proj_zarr_convention(obj: xr.Dataset | xr.DataArray) -> CRS | None:

tests/test_exact.py

@@ -6,6 +6,7 @@
 import hypothesis.strategies as st
 import numpy as np
 import numpy.testing as npt
+import pytest
 import sparse
 import xarray as xr
 import xarray.testing as xrt
@@ -142,3 +143,16 @@ def test_coverage_weights(
     if "area_spherical" not in coverage_weight:
         assert expected.nnz == actual_sparse.nnz
     xrt.assert_equal(dataset["spatial_ref"], actual["spatial_ref"])
+
+
+@pytest.mark.parametrize("chunked", [False, True])
+def test_coverage_ascending_y(chunked):
+    # dataset has descending y; flip it and expect the flipped result.
+    flipped = dataset.isel(y=slice(None, None, -1))
+    if chunked:
+        flipped = flipped.chunk(y=100, x=200)
+
+    expected = coverage(dataset, world[["geometry"]]).isel(y=slice(None, None, -1))
+    actual = coverage(flipped, world[["geometry"]])
+
+    npt.assert_equal(expected.data.todense(), actual.compute().data.todense())

tests/test_raster_index.py

@@ -834,10 +834,30 @@ def test_assign_index_from_coords():
     # Verify the transform
     # Coordinates are centered at pixels, so we expect identity transform at pixel corners
     # x[0] = 0.5, dx = 1.0 -> c = 0.5 - 1.0/2 = 0.0
-    # y[0] = 0.5, dy = 1.0 -> f = 0.5 - 1.0/2 = 0.0 but since y increases down, f = y[-1] - dy/2 = 9.5 - 0.5 = 9.0
-    expected_affine = Affine.translation(0.0, 9.0) * Affine.scale(1.0, 1.0)
+    # y[0] = 0.5, dy = 1.0 -> f = 0.5 - 1.0/2 = 0.0
+    expected_affine = Affine.translation(0.0, 0.0) * Affine.scale(1.0, 1.0)
     actual_affine = result.xindexes["x"].transform()
     assert actual_affine == expected_affine
+    xr.testing.assert_identical(result.x.variable.to_base_variable(), da.x.variable.to_base_variable())
+    xr.testing.assert_identical(result.y.variable.to_base_variable(), da.y.variable.to_base_variable())
+
+
+@pytest.mark.parametrize("y_ascending", [True, False])
+def test_assign_index_from_coords_roundtrip(y_ascending):
+    lon = np.arange(-179.5, 180, 90.0)
+    lat = np.array([-60.0, 0.0, 60.0])
+    if not y_ascending:
+        lat = lat[::-1]
+    ds = xr.Dataset({"foo": (("lat", "lon"), np.zeros((3, 4)))}, coords={"lat": lat, "lon": lon})
+
+    result = assign_index(ds, x_dim="lon", y_dim="lat")
+
+    np.testing.assert_array_equal(result.lat.values, lat)
+    np.testing.assert_array_equal(result.lon.values, lon)
+    affine = result.xindexes["lat"].transform()
+    assert affine == Affine.translation(-224.5, lat[0] - (lat[1] - lat[0]) / 2) * Affine.scale(
+        90.0, lat[1] - lat[0]
+    )
 
 
 def test_assign_index_dataset():

tests/test_rasterize.py

@@ -154,3 +154,19 @@ def test_geometry_mask_with_raster_index(engine, raster_index_dataset, dataset):
     xr.testing.assert_equal(result, expected)
     assert isinstance(result.xindexes["longitude"], RasterIndex)
     assert isinstance(result.xindexes["latitude"], RasterIndex)
+
+
+@pytest.mark.parametrize("chunked", [False, True])
+@pytest.mark.parametrize("func", [rasterize, geometry_mask])
+def test_ascending_y(engine, dataset, func, chunked):
+    # dataset has descending latitude; flip it and expect the flipped result.
+    world = gpd.read_file(geodatasets.get_path("naturalearth land"))
+    kwargs = dict(xdim="longitude", ydim="latitude", engine=engine)
+    flipped = dataset.isel(latitude=slice(None, None, -1))
+    if chunked:
+        flipped = flipped.chunk(latitude=100, longitude=200)
+
+    expected = func(dataset, world[["geometry"]], **kwargs).isel(latitude=slice(None, None, -1))
+    result = func(flipped, world[["geometry"]], **kwargs)
+
+    xr.testing.assert_equal(result.compute(), expected)