feat: Implementation of resample_spatial to reproject values

openeo_processes_dask/process_implementations/cubes/resample.py

@@ -0,0 +1,79 @@
+import logging
+from typing import Union
+
+import odc.geo
+from odc.geo.geobox import resolution_from_affine
+from pyproj.crs import CRS, CRSError
+
+from openeo_processes_dask.process_implementations.data_model import RasterCube
+
+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: Union[str, int] = None,
+    resolution: Union[str, int] = None,
+    method: str = "near",
+    align: str = "upper-left",
+):
+    """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."""
+
+    if align == "upper-left":
+        logging.warning("Warning: align parameter is unused by current implementation.")
+
+    # Assert resampling method is correct.
+    if method == "near":
+        method = "nearest"
+
+    elif method not in resample_methods_list:
+        raise Exception(
+            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("bands", "t", "y", "x")
+
+    if not projection:
+        projection = data_cp.rio.crs
+
+    try:
+        projection = CRS(projection)
+    except CRSError as e:
+        raise CRSError(f"{projection} Can not be parsed to CRS.")
+
+    if not resolution:
+        resolution = resolution_from_affine(data_cp.geobox.affine).x
+
+    reprojected = data_cp.odc.reproject(
+        how=projection, resolution=resolution, resampling=method
+    )
+
+    if "longitude" in reprojected.dims:
+        reprojected = reprojected.rename({"longitude": "x"})
+
+    if "latitude" in reprojected.dims:
+        reprojected = reprojected.rename({"latitude": "y"})
+
+    reprojected.attrs["crs"] = data_cp.rio.crs
+
+    # Undo odc specific re-ordering.
+    reprojected = reprojected.transpose(*data.dims)
+
+    return reprojected

openeo_processes_dask/process_implementations/cubes/utils.py

@@ -1,8 +1,8 @@
 try:
     import dask
-
 except ImportError:
     dask = None
+import math
 
 from xarray.core.duck_array_ops import isnull as xr_isnull
 

pyproject.toml

@@ -32,7 +32,7 @@ xgboost = { version = "^1.5.1", optional = true }
 rioxarray = { version = ">=0.12.0,<1", optional = true }
 odc-algo = { version = "==0.2.3", optional = true }
 openeo-pg-parser-networkx = { version = ">=2023.5.1", optional = true }
-odc-geo = { version = "^0.3.2", optional = true }
+odc-geo = { version = "^0.4.0", optional = true }
 
 [tool.poetry.group.dev.dependencies]
 pytest = "^7.2.0"

tests/test_resample.py

@@ -0,0 +1,60 @@
+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",
+    )
+
+    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.geobox.affine).x == output_res
+        assert resolution_from_affine(output_cube.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