From 022ebc04bef2cf6a03d0ba7fb31c623727725355 Mon Sep 17 00:00:00 2001 From: ValentinaHutter <85164505+ValentinaHutter@users.noreply.github.com> Date: Tue, 26 Nov 2024 13:58:48 +0100 Subject: [PATCH] Resample cube update (#300) * draft: test * fix resample cube spatial, not crop data * pre-commit run --- .../process_implementations/cubes/resample.py | 113 +++++++----------- tests/test_resample.py | 42 +++++++ 2 files changed, 88 insertions(+), 67 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index a9d80b0..c19f5a5 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -43,6 +43,30 @@ def resample_spatial( ): """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 data.openeo.y_dim is None or data.openeo.x_dim is None: + raise DimensionMissing(f"Spatial dimension missing for dataset: {data} ") + + methods_list = [ + "near", + "bilinear", + "cubic", + "cubicspline", + "lanczos", + "average", + "mode", + "max", + "min", + "med", + "q1", + "q3", + ] + + if method not in methods_list: + raise Exception( + f'Selected resampling method "{method}" is not available! Please select one of ' + f"[{', '.join(methods_list)}]" + ) + # Assert resampling method is correct. if method == "near": method = "nearest" @@ -53,18 +77,9 @@ def resample_spatial( f"[{', '.join(resample_methods_list)}]" ) - dims = list(data.dims) + dim_order = data.dims - known_dims = [] - if len(data.openeo.band_dims) > 0: - known_dims.append(data.openeo.band_dims[0]) - if len(data.openeo.temporal_dims) > 0: - known_dims.append(data.openeo.temporal_dims[0]) - known_dims += [data.openeo.y_dim, data.openeo.x_dim] - - other_dims = [dim for dim in dims if dim not in known_dims] - - data_cp = data.transpose(*other_dims, *known_dims) + data_cp = data.transpose(..., data.openeo.y_dim, data.openeo.x_dim) if projection is None: projection = data_cp.rio.crs @@ -89,6 +104,8 @@ def resample_spatial( if reprojected.openeo.y_dim != data.openeo.y_dim: reprojected = reprojected.rename({reprojected.openeo.y_dim: data.openeo.y_dim}) + reprojected = reprojected.transpose(*dim_order) + reprojected.attrs["crs"] = data_cp.rio.crs return reprojected @@ -97,67 +114,29 @@ def resample_spatial( def resample_cube_spatial( data: RasterCube, target: RasterCube, method="near", options=None ) -> RasterCube: - methods_list = [ - "near", - "bilinear", - "cubic", - "cubicspline", - "lanczos", - "average", - "mode", - "max", - "min", - "med", - "q1", - "q3", - ] - - if ( - data.openeo.y_dim is None - or data.openeo.x_dim is None - or target.openeo.y_dim is None - or target.openeo.x_dim is None - ): + if target.openeo.y_dim is None or target.openeo.x_dim is None: raise DimensionMissing( - f"Spatial dimension missing from data or target. Available dimensions for data: {data.dims} for target: {target.dims}" + f"Spatial dimension missing for target dataset: {target} " ) - # ODC reproject requires y to be before x - required_dim_order = (..., data.openeo.y_dim, data.openeo.x_dim) - - data_reordered = data.transpose(*required_dim_order, missing_dims="ignore") - target_reordered = target.transpose(*required_dim_order, missing_dims="ignore") - - if method == "near": - method = "nearest" - - elif method not in methods_list: - raise Exception( - f'Selected resampling method "{method}" is not available! Please select one of ' - f"[{', '.join(methods_list)}]" - ) - - resampled_data = data_reordered.odc.reproject( - target_reordered.odc.geobox, resampling=method + target_resolution, target_crs = None, None + if hasattr(target, "rio"): + if hasattr(target.rio, "resolution"): + if type(target.rio.resolution()) in [tuple, list]: + target_resolution = target.rio.resolution()[0] + else: + target_resolution = target.rio.resolution() + if hasattr(target.rio, "crs"): + target_crs = target.rio.crs + if not target_crs: + raise OpenEOException(f"Projection not found in target dataset: {target} ") + if not target_resolution: + raise OpenEOException(f"Resolution not found in target dataset: {target} ") + + resampled_data = resample_spatial( + data=data, projection=target_crs, resolution=target_resolution, method=method ) - resampled_data.rio.write_crs(target_reordered.rio.crs, inplace=True) - - try: - # odc.reproject renames the coordinates according to the geobox, this undoes that. - resampled_data = resampled_data.rename( - {"longitude": data.openeo.x_dim, "latitude": data.openeo.y_dim} - ) - except ValueError: - pass - - # Order axes back to how they were before - resampled_data = resampled_data.transpose(*data.dims) - - # Ensure that attrs except crs are copied over - for k, v in data.attrs.items(): - if k.lower() != "crs": - resampled_data.attrs[k] = v return resampled_data diff --git a/tests/test_resample.py b/tests/test_resample.py index 7efc7bb..0ef21a2 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -140,6 +140,48 @@ def test_resample_cube_spatial( assert output_cube.odc.spatial_dims == ("y", "x") +@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", + ], +) +@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_cube_spatial_small( + 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", + ) + + resampled_cube = resample_spatial( + data=input_cube, projection=output_crs, resolution=output_res + ) + + output_cube = resample_cube_spatial( + data=input_cube, target=resampled_cube[10:60, 20:150, :, :], method="average" + ) + + general_output_checks( + input_cube=input_cube, + output_cube=output_cube, + expected_dims=input_cube.dims, + verify_attrs=False, + verify_crs=False, + ) + + assert list(output_cube.shape) == list(resampled_cube.shape) + + @pytest.mark.parametrize("size", [(6, 5, 30, 4)]) @pytest.mark.parametrize("dtype", [np.float64]) @pytest.mark.parametrize(