resample_cube_spatial: align resampled coordinates (#306)

* add update to align resampled coordinates

* run pre-commit

* update resample

* check spatial coordinates in resample

* check spatial coordinates in resample

openeo_processes_dask/process_implementations/cubes/resample.py

@@ -137,6 +137,68 @@ def resample_cube_spatial(
         data=data, projection=target_crs, resolution=target_resolution, method=method
     )
 
+    def align_arrays(a, b, res):
+        desc = 1
+        if len(a) == 1 and len(b) == 1:
+            return a
+        elif len(a) == 1 and len(b) > 1:
+            if b[0] > b[1]:
+                desc = -1
+        elif len(a) > len(b):
+            b0 = np.absolute(b[0] - a).argmin()
+            blen = len(b)
+            close = np.isclose(b, a[b0 : b0 + blen], atol=res * 0.99)
+            if close.all():
+                return a[b0 : b0 + blen]
+            else:
+                raise Exception("Coordinates could not be aligned! ")
+        elif len(a) == len(b):
+            if b[0] > b[1]:
+                if a[0] < a[1]:
+                    a = np.flip(a)
+            else:
+                if a[0] > a[1]:
+                    a = np.flip(a)
+            close = np.isclose(b, a, atol=res * 0.99)
+            if close.all():
+                return a
+            else:
+                raise Exception("Coordinates could not be aligned! ")
+        else:
+            if b[0] > b[1]:
+                desc = -1
+                if a[0] < a[1]:
+                    a = np.flip(a)
+            else:
+                if a[0] > a[1]:
+                    a = np.flip(a)
+        a0 = np.absolute(b - a[0]).argmin()
+        alen = len(a)
+        close = np.isclose(a, b[a0 : a0 + alen], atol=res * 0.99)
+        if not close.all():
+            raise Exception("Coordinates could not be aligned! ")
+
+        new_b = np.arange(b[0], a[0], res * desc)
+        new_b = np.append(new_b, a)
+        new_b = np.append(new_b, np.flip(np.arange(b[-1], a[-1], -res * desc)))
+        if len(b) != len(new_b):
+            raise Exception("Coordinates could not be aligned! ")
+        return new_b
+
+    x_data = resampled_data[resampled_data.openeo.x_dim[0]].values
+    x_target = target[target.openeo.x_dim[0]].values
+    if not (len(x_data) == len(x_target) and (x_data == x_target).all()):
+        resampled_data[resampled_data.openeo.x_dim[0]] = align_arrays(
+            x_target, x_data, target_resolution
+        )
+
+    y_data = resampled_data[resampled_data.openeo.y_dim[0]].values
+    y_target = target[target.openeo.y_dim[0]].values
+    if not (len(y_data) == len(y_target) and (y_data == y_target).all()):
+        resampled_data[resampled_data.openeo.y_dim[0]] = align_arrays(
+            y_target, y_data, target_resolution
+        )
+
     return resampled_data
 
 

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "openeo-processes-dask"
-version = "2024.11.5"
+version = "2024.12.1"
 description = "Python implementations of many OpenEO processes, dask-friendly by default."
 authors = ["Lukas Weidenholzer <[email protected]>", "Sean Hoyal <[email protected]>", "Valentina Hutter <[email protected]>"]
 maintainers = ["EODC Staff <[email protected]>"]

tests/test_resample.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+import xarray as xr
 from odc.geo.geobox import resolution_from_affine
 from openeo_pg_parser_networkx.pg_schema import ParameterReference, TemporalInterval
 from pyproj.crs import CRS
@@ -180,6 +181,8 @@ def test_resample_cube_spatial_small(
     )
 
     assert list(output_cube.shape) == list(resampled_cube.shape)
+    assert (output_cube["x"].values == resampled_cube["x"].values).all()
+    assert (output_cube["y"].values == resampled_cube["y"].values).all()
 
 
 @pytest.mark.parametrize("size", [(6, 5, 30, 4)])