feat: Implementation of resample_spatial to reproject values (#103)

* basic implementation of resample_spatial to reproject values

* update function to handle spatial resampling as well as reprojection

* fix from pre commit

* pre commit fix

* pre commit file edits

* accept reformatting

* add reformatting

* add reformat

* accept reformats

* accept reformatting

* drop broken import

* accept reformatting

* accept reformatting

* fixed broken tests

* add comment

* add negative test case

* resolve locally

* try merge from main to fix tests and change crs function

* Update openeo_processes_dask/process_implementations/cubes/resample.py

Co-authored-by: Lukas Weidenholzer <[email protected]>

* Update openeo_processes_dask/process_implementations/cubes/resample.py

Co-authored-by: Lukas Weidenholzer <[email protected]>

* Update openeo_processes_dask/process_implementations/cubes/resample.py

Co-authored-by: Lukas Weidenholzer <[email protected]>

* Update openeo_processes_dask/process_implementations/cubes/resample.py

Co-authored-by: Lukas Weidenholzer <[email protected]>

* comments updates

* import typing optional

* fix tessts

* bump process specs

---------

Co-authored-by: Lukas Weidenholzer <[email protected]>
Co-authored-by: Lukas Weidenholzer <[email protected]>

openeo_processes_dask/process_implementations/cubes/resample.py

@@ -0,0 +1,80 @@
+import logging
+from typing import Optional, Union
+
+import odc.geo.xr
+from odc.geo.geobox import resolution_from_affine
+from pyproj.crs import CRS, CRSError
+
+from openeo_processes_dask.process_implementations.data_model import RasterCube
+from openeo_processes_dask.process_implementations.exceptions import OpenEOException
+
+logger = logging.getLogger(__name__)
+
+resample_methods_list = [
+    "near",
+    "bilinear",
+    "cubic",
+    "cubicspline",
+    "lanczos",
+    "average",
+    "mode",
+    "max",
+    "min",
+    "med",
+    "q1",
+    "q3",
+]
+
+
+def resample_spatial(
+    data: RasterCube,
+    projection: Optional[Union[str, int]] = None,
+    resolution: int = 0,
+    method: str = "near",
+):
+    """Resamples the spatial dimensions (x,y) of the data cube to a specified resolution and/or warps the data cube to the target projection. At least resolution or projection must be specified."""
+
+    # Assert resampling method is correct.
+    if method == "near":
+        method = "nearest"
+
+    elif method not in resample_methods_list:
+        raise OpenEOException(
+            f'Selected resampling method "{method}" is not available! Please select one of '
+            f"[{', '.join(resample_methods_list)}]"
+        )
+
+    # Re-order, this is specifically done for odc reproject
+    data_cp = data.transpose(
+        data.openeo.band_dims[0],
+        data.openeo.temporal_dims[0],
+        data.openeo.y_dim,
+        data.openeo.x_dim,
+    )
+
+    if projection is None:
+        projection = data_cp.rio.crs
+
+    try:
+        projection = CRS.from_user_input(projection)
+    except CRSError as e:
+        raise CRSError(
+            f"Provided projection string: '{projection}' can not be parsed to CRS."
+        ) from e
+
+    if resolution == 0:
+        resolution = resolution_from_affine(data_cp.odc.geobox.affine).x
+
+    reprojected = data_cp.odc.reproject(
+        how=projection, resolution=resolution, resampling=method
+    )
+
+    if "longitude" in reprojected.dims and "x" in data.dims:
+        reprojected = reprojected.rename({"longitude": "x"})
+
+    if "latitude" in reprojected.dims and "y" in data.dims:
+        reprojected = reprojected.rename({"latitude": "y"})
+
+    reprojected.attrs["crs"] = data_cp.rio.crs
+
+    return reprojected

openeo_processes_dask/process_implementations/cubes/utils.py

@@ -1,6 +1,5 @@
 try:
     import dask
-
 except ImportError:
     dask = None
 

pyproject.toml

@@ -32,7 +32,7 @@ dask-geopandas = { version = ">=0.2.0,<1", optional = true }
 xgboost = { version = "^1.5.1", optional = true }
 rioxarray = { version = ">=0.12.0,<1", optional = true }
 openeo-pg-parser-networkx = { version = ">=2023.5.1", optional = true }
-odc-geo = { version = ">=0.3.2,<1", optional = true }
+odc-geo = { version = ">=0.4.1,<1", optional = true }
 stac_validator = { version = ">=3.3.1", optional = true }
 stackstac = { version = ">=0.4.3", optional = true }
 pystac_client = { version = ">=0.6.1", optional = true }

tests/test_resample.py

@@ -0,0 +1,65 @@
+import numpy as np
+import pytest
+from odc.geo.geobox import resolution_from_affine
+from pyproj.crs import CRS
+
+from openeo_processes_dask.process_implementations.cubes.resample import (
+    resample_spatial,
+)
+from tests.general_checks import general_output_checks
+from tests.mockdata import create_fake_rastercube
+
+
+@pytest.mark.parametrize(
+    "output_crs",
+    [
+        3587,
+        "32633",
+        "+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.81977 +y_0=2121415.69617 +datum=WGS84 +units=m +no_defs",
+        "4326",
+    ],
+)
+@pytest.mark.parametrize("output_res", [5, 30, 60])
+@pytest.mark.parametrize("size", [(30, 30, 20, 4)])
+@pytest.mark.parametrize("dtype", [np.float32])
+def test_resample_spatial(
+    output_crs, output_res, temporal_interval, bounding_box, random_raster_data
+):
+    """Test to ensure resolution gets changed correctly."""
+    input_cube = create_fake_rastercube(
+        data=random_raster_data,
+        spatial_extent=bounding_box,
+        temporal_extent=temporal_interval,
+        bands=["B02", "B03", "B04", "B08"],
+        backend="dask",
+    )
+
+    with pytest.raises(Exception):
+        output_cube = resample_spatial(
+            data=input_cube, projection=output_crs, resolution=output_res, method="bad"
+        )
+
+    output_cube = resample_spatial(
+        data=input_cube, projection=output_crs, resolution=output_res
+    )
+
+    general_output_checks(
+        input_cube=input_cube,
+        output_cube=output_cube,
+        verify_attrs=False,
+        verify_crs=False,
+    )
+
+    assert output_cube.odc.spatial_dims == ("y", "x")
+    assert output_cube.rio.crs == CRS.from_user_input(output_crs)
+
+    if output_crs != "4326":
+        assert resolution_from_affine(output_cube.odc.geobox.affine).x == output_res
+        assert resolution_from_affine(output_cube.odc.geobox.affine).y == -output_res
+
+    if output_cube.rio.crs.is_geographic:
+        assert min(output_cube.x) >= -180
+        assert max(output_cube.x) <= 180
+
+        assert min(output_cube.y) >= -90
+        assert max(output_cube.y) <= 90