diff --git a/src/rasterix/rasterize/core.py b/src/rasterix/rasterize/core.py index 96aa2c8..6ca5930 100644 --- a/src/rasterix/rasterize/core.py +++ b/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) diff --git a/src/rasterix/rasterize/exact.py b/src/rasterix/rasterize/exact.py index e363413..6ae6e38 100644 --- a/src/rasterix/rasterize/exact.py +++ b/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) # ============================================================================ diff --git a/src/rasterix/utils.py b/src/rasterix/utils.py index 8a080b3..f3fce6e 100644 --- a/src/rasterix/utils.py +++ b/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: diff --git a/tests/test_exact.py b/tests/test_exact.py index 5b8fb20..81ff943 100644 --- a/tests/test_exact.py +++ b/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()) diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py index a77e973..32f3d0e 100644 --- a/tests/test_raster_index.py +++ b/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(): diff --git a/tests/test_rasterize.py b/tests/test_rasterize.py index e826d05..8af6050 100644 --- a/tests/test_rasterize.py +++ b/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)