From d1861343bcc2a137e0e471179c7148deb7a1147f Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 4 May 2023 16:26:54 +0200 Subject: [PATCH 01/26] basic implementation of resample_spatial to reproject values --- .../process_implementations/cubes/resample.py | 94 +++++++++++++++++++ .../process_implementations/cubes/utils.py | 34 ++++++- tests/test_resample.py | 34 +++++++ 3 files changed, 160 insertions(+), 2 deletions(-) create mode 100644 openeo_processes_dask/process_implementations/cubes/resample.py create mode 100644 tests/test_resample.py diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py new file mode 100644 index 00000000..bab4e9b2 --- /dev/null +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -0,0 +1,94 @@ +import odc.algo + +from affine import Affine +from enum import Enum +from datacube.utils.geometry import GeoBox +from math import pi +from pyproj import Transformer +from pyproj.crs import CRS, CRSError + +from typing import Union, Optional + +from openeo_processes_dask.process_implementations.cubes.utils import detect_changing_unit +from openeo_processes_dask.process_implementations.data_model import RasterCube + +__all__ = ["resample_spatial"] + + +def resample_spatial( + data: RasterCube, + epsg_code: Union[str, int], + resolution: int = 0, + method: str = 'near', + align: str = "upper-left"): + """ Resample the input rastercube by the provided epsg_code. Currently not accepting resolution resampling. """ + # Re-order + data = data.transpose("bands", "t", "y", "x") + + src_crs = CRS(data.rio.crs) + + try: + dst_crs = CRS.from_epsg(epsg_code) + except CRSError: + raise Exception(f"epsg_code parameter {epsg_code} is not a valid epsg code.") + + transformer = Transformer.from_crs( + src_crs, + dst_crs, + always_xy=True, + ) + + # For the affine, we want to know the x and y coordinate for the upper-left pixel. + tmp_left_pixel = data.isel( + x=[0], + y=[0] + ) + tmp_left_pixel.x.values, tmp_left_pixel.y.values + # Transform the upper left pixel coords to the destination crs. + new_x, new_y = transformer.transform( + tmp_left_pixel.x.values, + tmp_left_pixel.y.values + ) + + resolution = detect_changing_unit( + src_crs=src_crs, + dst_crs= dst_crs, + # First value of affine is pixel width. + src_res=data.affine[0] + ) + + # Docs for geotransform + # https://gdal.org/tutorials/geotransforms_tut.html + new_affine = Affine( + resolution, + 0, + new_x[0], + 0, + # Negative pixel width used for pixel height + -resolution, + new_y[0] + ) + + dst_geo = GeoBox( + width=len(data.x), + height=len(data.y), + crs=dst_crs, + affine=new_affine + ) + + reprojected = odc.algo.xr_reproject( + src=data, + geobox=dst_geo, + ) + if "longitude" in reprojected.dims: + reprojected = reprojected.rename( + {"longitude": data.openeo.x_dim} + ) + + if "latitude" in reprojected.dims: + reprojected = reprojected.rename( + {"latitude": data.openeo.y_dim} + ) + + reprojected.attrs["crs"] = dst_crs + return reprojected \ No newline at end of file diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index fa5a06e1..868da580 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -1,9 +1,9 @@ try: import dask - except ImportError: dask = None - +from enum import Enum +from math import pi def _has_dask(): return dask is not None @@ -11,3 +11,33 @@ def _has_dask(): def _is_dask_array(arr): return _has_dask() and isinstance(arr, dask.array.Array) + + +class Unit(Enum): + metre = 'metre' + degree = 'degree' + +def approx_metres_2_degrees(length_metres = 10): + """ Angle of rotation calculated off the earths equatorial radius in metres. + Ref: https://en.ans.wiki/5729/how-to-convert-meters-to-degrees/ """ + earth_equator_radius_metres = 6378160 + return (180 * length_metres ) / (earth_equator_radius_metres * pi) + +def approx_degrees_2_metres(angle_in_degrees = 0.0009): + """ Calculated the metres travelled for an angle of rotation on the earths equatorial radius in metres. + Ref: https://en.ans.wiki/5728/how-to-convert-degrees-to-meters/ """ + earth_equator_radius_metres = 6378160 + return (earth_equator_radius_metres * pi * angle_in_degrees) / 180 + +def detect_changing_unit(src_crs, dst_crs, src_res): + """ Is there a unit change between the src and dst CRS. If so, return the approximate value for the new resolution. """ + + src_unit = Unit(src_crs.axis_info[0].unit_name) + dst_unit = Unit(dst_crs.axis_info[0].unit_name) + + if src_unit != dst_unit: + if dst_unit == Unit.metre: + return approx_degrees_2_metres(src_res) + if dst_unit == Unit.degree: + return approx_metres_2_degrees(src_res) + return src_res \ No newline at end of file diff --git a/tests/test_resample.py b/tests/test_resample.py new file mode 100644 index 00000000..ed7b8cf6 --- /dev/null +++ b/tests/test_resample.py @@ -0,0 +1,34 @@ +import numpy as np +import pytest + +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("size", [(30, 30, 20, 4)]) +@pytest.mark.parametrize("dtype", [np.float32]) +def test_resample_spatial( + temporal_interval, bounding_box, random_raster_data +): + 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, + epsg_code=32633 + ) + + general_output_checks( + input_cube=input_cube, + output_cube=output_cube, + verify_attrs=False, + verify_crs=False, + ) + + assert output_cube.rio.crs == CRS.from_epsg(32633) From 8d474fbba5f5bfe16940c9a957ff74ef8401944b Mon Sep 17 00:00:00 2001 From: sean Date: Fri, 12 May 2023 16:19:16 +0200 Subject: [PATCH 02/26] update function to handle spatial resampling as well as reprojection --- .../process_implementations/cubes/resample.py | 185 ++++++++++-------- .../process_implementations/cubes/utils.py | 71 +++++-- 2 files changed, 162 insertions(+), 94 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index bab4e9b2..f2466009 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -1,94 +1,119 @@ -import odc.algo +from typing import Union -from affine import Affine -from enum import Enum -from datacube.utils.geometry import GeoBox -from math import pi +import odc.algo from pyproj import Transformer from pyproj.crs import CRS, CRSError -from typing import Union, Optional - -from openeo_processes_dask.process_implementations.cubes.utils import detect_changing_unit +from openeo_processes_dask.process_implementations.cubes.utils import ( + detect_changing_unit, +) from openeo_processes_dask.process_implementations.data_model import RasterCube -__all__ = ["resample_spatial"] +resample_methods_list = [ + "near", + "bilinear", + "cubic", + "cubicspline", + "lanczos", + "average", + "mode", + "max", + "min", + "med", + "q1", + "q3", +] def resample_spatial( - data: RasterCube, - epsg_code: Union[str, int], - resolution: int = 0, - method: str = 'near', - align: str = "upper-left"): - """ Resample the input rastercube by the provided epsg_code. Currently not accepting resolution resampling. """ - # Re-order + data: RasterCube, + epsg_code: Union[str, int] = None, + resolution: int = None, + method: str = "near", + align: str = "upper-left", +): + """Resample the input rastercube by the provided epsg_code. Currently not accepting resolution resampling.""" + + # 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 xr_reproject data = data.transpose("bands", "t", "y", "x") - src_crs = CRS(data.rio.crs) - - try: - dst_crs = CRS.from_epsg(epsg_code) - except CRSError: - raise Exception(f"epsg_code parameter {epsg_code} is not a valid epsg code.") - - transformer = Transformer.from_crs( - src_crs, - dst_crs, - always_xy=True, - ) - - # For the affine, we want to know the x and y coordinate for the upper-left pixel. - tmp_left_pixel = data.isel( - x=[0], - y=[0] - ) - tmp_left_pixel.x.values, tmp_left_pixel.y.values - # Transform the upper left pixel coords to the destination crs. - new_x, new_y = transformer.transform( - tmp_left_pixel.x.values, - tmp_left_pixel.y.values - ) - - resolution = detect_changing_unit( - src_crs=src_crs, - dst_crs= dst_crs, - # First value of affine is pixel width. - src_res=data.affine[0] - ) - - # Docs for geotransform - # https://gdal.org/tutorials/geotransforms_tut.html - new_affine = Affine( - resolution, - 0, - new_x[0], - 0, - # Negative pixel width used for pixel height - -resolution, - new_y[0] - ) - - dst_geo = GeoBox( - width=len(data.x), - height=len(data.y), - crs=dst_crs, - affine=new_affine - ) - - reprojected = odc.algo.xr_reproject( - src=data, - geobox=dst_geo, - ) - if "longitude" in reprojected.dims: - reprojected = reprojected.rename( - {"longitude": data.openeo.x_dim} + # Do reprojection first, and then resampling + if epsg_code: + try: + dst_crs = CRS.from_epsg(epsg_code) + except CRSError: + raise Exception( + f"epsg_code parameter {epsg_code} is not a valid epsg code." + ) + + # Get original crs and resolution from dataset. + src_crs = CRS(data.rio.crs) + + transformer = Transformer.from_crs( + src_crs, + dst_crs, + always_xy=True, + ) + # For the affine, we want to know the x and y coordinate for the upper-left pixel. + top_left_pixel = data.isel(x=[0], y=[0]) + + # Transform the upper left pixel coords to the destination crs. + new_x, new_y = transformer.transform( + top_left_pixel.x.values, top_left_pixel.y.values ) + # If resolution has not been provided, it must be inferred from old resolution. + # First value of affine is pixel width. Use it as "resolution". + src_res = data.affine[0] + dst_res = detect_changing_unit( + src_crs=src_crs, + dst_crs=dst_crs, + # First value of affine is pixel width, which we use as 'resolution'. + src_res=src_res, + ) + + dst_geobox = prepare_geobox(data, dst_crs, dst_res, src_res, new_x[0], new_y[0]) + data = odc.algo.xr_reproject(src=data, geobox=dst_geobox) + data.attrs["crs"] = dst_crs + + # And if resampling was requested + if resolution: + dst_res = resolution + # Get src_crs seperately incase reprojection was carried out. + src_crs = CRS(data.rio.crs) + # Get top left pixel seperately incase reprojection was carried out + top_left_pixel = data.isel(x=[0], y=[0]) + src_res = data.affine[0] + + dst_geobox = prepare_geobox( + data, + src_crs, + dst_res, + src_res, + top_left_pixel.x.values[0], + top_left_pixel.y.values[0], + scale=True, + ) + + data = odc.algo.xr_reproject(src=data, geobox=dst_geobox) + + reprojected = data + + if "longitude" in reprojected.dims: + reprojected = reprojected.rename({"longitude": data.openeo.x_dim}) + if "latitude" in reprojected.dims: - reprojected = reprojected.rename( - {"latitude": data.openeo.y_dim} - ) - - reprojected.attrs["crs"] = dst_crs - return reprojected \ No newline at end of file + reprojected = reprojected.rename({"latitude": data.openeo.y_dim}) + + reprojected.attrs["crs"] = data.rio.crs + return reprojected diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index 868da580..e1a98cb4 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -2,8 +2,12 @@ import dask except ImportError: dask = None +import math from enum import Enum -from math import pi + +from affine import Affine +from odc.geo.geobox import GeoBox, resolution_from_affine + def _has_dask(): return dask is not None @@ -14,30 +18,69 @@ def _is_dask_array(arr): class Unit(Enum): - metre = 'metre' - degree = 'degree' + metre = "metre" + degree = "degree" -def approx_metres_2_degrees(length_metres = 10): - """ Angle of rotation calculated off the earths equatorial radius in metres. - Ref: https://en.ans.wiki/5729/how-to-convert-meters-to-degrees/ """ + +def approx_metres_2_degrees(length_metres=10): + """Angle of rotation calculated off the earths equatorial radius in metres. + Ref: https://en.ans.wiki/5729/how-to-convert-meters-to-degrees/""" earth_equator_radius_metres = 6378160 - return (180 * length_metres ) / (earth_equator_radius_metres * pi) + return ((180 * length_metres) / (earth_equator_radius_metres * math.pi)) * 10 + -def approx_degrees_2_metres(angle_in_degrees = 0.0009): - """ Calculated the metres travelled for an angle of rotation on the earths equatorial radius in metres. - Ref: https://en.ans.wiki/5728/how-to-convert-degrees-to-meters/ """ +def approx_degrees_2_metres(angle_in_degrees=0.0009): + """Calculated the metres travelled for an angle of rotation on the earths equatorial radius in metres. + Ref: https://en.ans.wiki/5728/how-to-convert-degrees-to-meters/""" earth_equator_radius_metres = 6378160 - return (earth_equator_radius_metres * pi * angle_in_degrees) / 180 + return ((earth_equator_radius_metres * math.pi * angle_in_degrees) / 180) / 10 -def detect_changing_unit(src_crs, dst_crs, src_res): - """ Is there a unit change between the src and dst CRS. If so, return the approximate value for the new resolution. """ + +def detect_changing_unit(src_crs, dst_crs, src_res=None): + """Is there a unit change between the src and dst CRS. If so, return the approximate value for the new resolution.""" src_unit = Unit(src_crs.axis_info[0].unit_name) dst_unit = Unit(dst_crs.axis_info[0].unit_name) + # If not dst_res has been set, check the src_unit and dst_unit to + # determine whether the src_res needs to be converted. if src_unit != dst_unit: if dst_unit == Unit.metre: + print("deg 2 met") return approx_degrees_2_metres(src_res) if dst_unit == Unit.degree: + print("met 2 deg") return approx_metres_2_degrees(src_res) - return src_res \ No newline at end of file + return src_res + + +def prepare_geobox(data, dst_crs, dst_res, new_top_left_x, new_top_left_y, scale=False): + """Get the destination geobox ready for the transformation.""" + # Docs for geotransform + # https://gdal.org/tutorials/geotransforms_tut.html + + new_affine = Affine( + dst_res, + 0, + new_top_left_x, + 0, + # Negative pixel width used for pixel height + -dst_res, + new_top_left_y, + ) + + if scale: + updated_resolution = resolution_from_affine(new_affine) + + x_min, x_max = min(data.x.values), max(data.x.values) + y_min, y_max = min(data.y.values), max(data.y.values) + + x_length = math.ceil(abs(x_min - x_max) / updated_resolution.x) + y_length = math.ceil(abs(y_min - y_max) / -updated_resolution.y) + + else: + x_length = len(data.x) + y_length = len(data.y) + + dst_geo = GeoBox(width=x_length, height=y_length, crs=dst_crs, affine=new_affine) + return dst_geo From 87680cc2a4829a420a872a9a1d5d625180e2d550 Mon Sep 17 00:00:00 2001 From: sean Date: Mon, 15 May 2023 08:49:30 +0200 Subject: [PATCH 03/26] fix from pre commit --- tests/test_resample.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/tests/test_resample.py b/tests/test_resample.py index ed7b8cf6..e369360d 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -1,16 +1,17 @@ import numpy as np import pytest - from pyproj.crs import CRS -from openeo_processes_dask.process_implementations.cubes.resample import resample_spatial + +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("size", [(30, 30, 20, 4)]) @pytest.mark.parametrize("dtype", [np.float32]) -def test_resample_spatial( - temporal_interval, bounding_box, random_raster_data -): +def test_resample_spatial(temporal_interval, bounding_box, random_raster_data): input_cube = create_fake_rastercube( data=random_raster_data, spatial_extent=bounding_box, @@ -19,10 +20,7 @@ def test_resample_spatial( backend="dask", ) - output_cube = resample_spatial( - data=input_cube, - epsg_code=32633 - ) + output_cube = resample_spatial(data=input_cube, epsg_code=32633) general_output_checks( input_cube=input_cube, From c4c77237e47978497850529b46b49baa623bb6a2 Mon Sep 17 00:00:00 2001 From: sean Date: Mon, 15 May 2023 09:20:46 +0200 Subject: [PATCH 04/26] pre commit fix --- .../process_implementations/cubes/resample.py | 1 + openeo_processes_dask/process_implementations/cubes/utils.py | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index f2466009..46c2c9e5 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -6,6 +6,7 @@ from openeo_processes_dask.process_implementations.cubes.utils import ( detect_changing_unit, + prepare_geobox, ) from openeo_processes_dask.process_implementations.data_model import RasterCube diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index e1a98cb4..085f3e77 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -6,7 +6,8 @@ from enum import Enum from affine import Affine -from odc.geo.geobox import GeoBox, resolution_from_affine +from datacube.utils.geometry import GeoBox +from odc.geo.geobox import resolution_from_affine def _has_dask(): From 48931e90f7bccb8da9bda59dbab590f089691ccb Mon Sep 17 00:00:00 2001 From: sean Date: Wed, 24 May 2023 10:06:28 +0200 Subject: [PATCH 05/26] pre commit file edits --- .../process_implementations/cubes/resample.py | 5 +++-- .../process_implementations/cubes/utils.py | 4 +++- tests/test_resample.py | 11 ++++++++--- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 46c2c9e5..65cbd41b 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -1,6 +1,7 @@ from typing import Union import odc.algo +from odc.geo.geobox import resolution_from_affine from pyproj import Transformer from pyproj.crs import CRS, CRSError @@ -75,7 +76,7 @@ def resample_spatial( # If resolution has not been provided, it must be inferred from old resolution. # First value of affine is pixel width. Use it as "resolution". - src_res = data.affine[0] + src_res = resolution_from_affine(data.affine).x dst_res = detect_changing_unit( src_crs=src_crs, dst_crs=dst_crs, @@ -94,7 +95,7 @@ def resample_spatial( src_crs = CRS(data.rio.crs) # Get top left pixel seperately incase reprojection was carried out top_left_pixel = data.isel(x=[0], y=[0]) - src_res = data.affine[0] + src_res = resolution_from_affine(data.affine).x dst_geobox = prepare_geobox( data, diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index 085f3e77..c35c89e7 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -55,7 +55,9 @@ def detect_changing_unit(src_crs, dst_crs, src_res=None): return src_res -def prepare_geobox(data, dst_crs, dst_res, new_top_left_x, new_top_left_y, scale=False): +def prepare_geobox( + data, dst_crs, dst_res, src_res, new_top_left_x, new_top_left_y, scale=False +): """Get the destination geobox ready for the transformation.""" # Docs for geotransform # https://gdal.org/tutorials/geotransforms_tut.html diff --git a/tests/test_resample.py b/tests/test_resample.py index e369360d..5a5469e4 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -9,9 +9,13 @@ from tests.mockdata import create_fake_rastercube +@pytest.mark.parametrize("output_crs", [3587, 32633, 6068]) @pytest.mark.parametrize("size", [(30, 30, 20, 4)]) @pytest.mark.parametrize("dtype", [np.float32]) -def test_resample_spatial(temporal_interval, bounding_box, random_raster_data): +def test_resample_spatial( + output_crs, temporal_interval, bounding_box, random_raster_data +): + """Test to ensure CRS get changed correctly.""" input_cube = create_fake_rastercube( data=random_raster_data, spatial_extent=bounding_box, @@ -20,7 +24,7 @@ def test_resample_spatial(temporal_interval, bounding_box, random_raster_data): backend="dask", ) - output_cube = resample_spatial(data=input_cube, epsg_code=32633) + output_cube = resample_spatial(data=input_cube, epsg_code=output_crs) general_output_checks( input_cube=input_cube, @@ -29,4 +33,5 @@ def test_resample_spatial(temporal_interval, bounding_box, random_raster_data): verify_crs=False, ) - assert output_cube.rio.crs == CRS.from_epsg(32633) + assert output_cube.rio.crs == CRS.from_epsg(output_crs) + assert len(output_cube.x) == len(input_cube.x) From dde83c085b59fb87339b9183b27c994d6d853ced Mon Sep 17 00:00:00 2001 From: sean Date: Wed, 24 May 2023 13:53:20 +0200 Subject: [PATCH 06/26] accept reformatting --- .../process_implementations/cubes/resample.py | 16 ++++++------- tests/test_resample.py | 24 ++++++++++++++----- 2 files changed, 26 insertions(+), 14 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 65cbd41b..ab012d7b 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -29,7 +29,7 @@ def resample_spatial( data: RasterCube, - epsg_code: Union[str, int] = None, + projection: Union[str, int] = None, resolution: int = None, method: str = "near", align: str = "upper-left", @@ -50,12 +50,12 @@ def resample_spatial( data = data.transpose("bands", "t", "y", "x") # Do reprojection first, and then resampling - if epsg_code: + if projection: try: - dst_crs = CRS.from_epsg(epsg_code) + dst_crs = CRS.from_user_input(projection) except CRSError: raise Exception( - f"epsg_code parameter {epsg_code} is not a valid epsg code." + f"Projection parameter {projection} could not be parsed to CRS." ) # Get original crs and resolution from dataset. @@ -93,6 +93,7 @@ def resample_spatial( dst_res = resolution # Get src_crs seperately incase reprojection was carried out. src_crs = CRS(data.rio.crs) + # Get top left pixel seperately incase reprojection was carried out top_left_pixel = data.isel(x=[0], y=[0]) src_res = resolution_from_affine(data.affine).x @@ -107,15 +108,14 @@ def resample_spatial( scale=True, ) - data = odc.algo.xr_reproject(src=data, geobox=dst_geobox) + data = odc.algo.xr_reproject(src=data, geobox=dst_geobox, resampling=method) reprojected = data - if "longitude" in reprojected.dims: - reprojected = reprojected.rename({"longitude": data.openeo.x_dim}) + reprojected = reprojected.rename({"longitude": "x"}) if "latitude" in reprojected.dims: - reprojected = reprojected.rename({"latitude": data.openeo.y_dim}) + reprojected = reprojected.rename({"latitude": "y"}) reprojected.attrs["crs"] = data.rio.crs return reprojected diff --git a/tests/test_resample.py b/tests/test_resample.py index 5a5469e4..5ec4258d 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -1,5 +1,6 @@ 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 ( @@ -9,13 +10,21 @@ from tests.mockdata import create_fake_rastercube -@pytest.mark.parametrize("output_crs", [3587, 32633, 6068]) +@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_spatial( - output_crs, temporal_interval, bounding_box, random_raster_data + output_crs, output_res, temporal_interval, bounding_box, random_raster_data ): - """Test to ensure CRS get changed correctly.""" + """Test to ensure resolution gets changed correctly.""" input_cube = create_fake_rastercube( data=random_raster_data, spatial_extent=bounding_box, @@ -24,7 +33,9 @@ def test_resample_spatial( backend="dask", ) - output_cube = resample_spatial(data=input_cube, epsg_code=output_crs) + output_cube = resample_spatial( + data=input_cube, projection=output_crs, resolution=output_res + ) general_output_checks( input_cube=input_cube, @@ -33,5 +44,6 @@ def test_resample_spatial( verify_crs=False, ) - assert output_cube.rio.crs == CRS.from_epsg(output_crs) - assert len(output_cube.x) == len(input_cube.x) + assert output_cube.rio.crs == CRS.from_user_input(output_crs) + assert resolution_from_affine(output_cube.geobox.affine).x == output_res + assert resolution_from_affine(output_cube.geobox.affine).y == -output_res From 92562d5686c1a179169248c1e2c3462f34ffa82b Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 25 May 2023 09:01:22 +0200 Subject: [PATCH 07/26] add reformatting --- .../process_implementations/cubes/resample.py | 78 ++++--------------- .../process_implementations/cubes/utils.py | 71 ----------------- pyproject.toml | 2 +- 3 files changed, 14 insertions(+), 137 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index ab012d7b..3a28e9d3 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -1,14 +1,9 @@ from typing import Union -import odc.algo +import odc.geo from odc.geo.geobox import resolution_from_affine -from pyproj import Transformer from pyproj.crs import CRS, CRSError -from openeo_processes_dask.process_implementations.cubes.utils import ( - detect_changing_unit, - prepare_geobox, -) from openeo_processes_dask.process_implementations.data_model import RasterCube resample_methods_list = [ @@ -30,7 +25,7 @@ def resample_spatial( data: RasterCube, projection: Union[str, int] = None, - resolution: int = None, + resolution: Union[str, int] = None, method: str = "near", align: str = "upper-left", ): @@ -49,68 +44,21 @@ def resample_spatial( # Re-order, this is specifically done for xr_reproject data = data.transpose("bands", "t", "y", "x") - # Do reprojection first, and then resampling - if projection: - try: - dst_crs = CRS.from_user_input(projection) - except CRSError: - raise Exception( - f"Projection parameter {projection} could not be parsed to CRS." - ) + if not projection: + projection = data.rio.crs - # Get original crs and resolution from dataset. - src_crs = CRS(data.rio.crs) + try: + projection = CRS(projection) + except CRSError as e: + raise CRSError(f"{projection} Can not be parsed to CRS.") - transformer = Transformer.from_crs( - src_crs, - dst_crs, - always_xy=True, - ) - # For the affine, we want to know the x and y coordinate for the upper-left pixel. - top_left_pixel = data.isel(x=[0], y=[0]) - - # Transform the upper left pixel coords to the destination crs. - new_x, new_y = transformer.transform( - top_left_pixel.x.values, top_left_pixel.y.values - ) - - # If resolution has not been provided, it must be inferred from old resolution. - # First value of affine is pixel width. Use it as "resolution". - src_res = resolution_from_affine(data.affine).x - dst_res = detect_changing_unit( - src_crs=src_crs, - dst_crs=dst_crs, - # First value of affine is pixel width, which we use as 'resolution'. - src_res=src_res, - ) - - dst_geobox = prepare_geobox(data, dst_crs, dst_res, src_res, new_x[0], new_y[0]) - data = odc.algo.xr_reproject(src=data, geobox=dst_geobox) - data.attrs["crs"] = dst_crs - - # And if resampling was requested - if resolution: - dst_res = resolution - # Get src_crs seperately incase reprojection was carried out. - src_crs = CRS(data.rio.crs) - - # Get top left pixel seperately incase reprojection was carried out - top_left_pixel = data.isel(x=[0], y=[0]) - src_res = resolution_from_affine(data.affine).x - - dst_geobox = prepare_geobox( - data, - src_crs, - dst_res, - src_res, - top_left_pixel.x.values[0], - top_left_pixel.y.values[0], - scale=True, - ) + if not resolution: + resolution = resolution_from_affine(data.geobox.affine).x - data = odc.algo.xr_reproject(src=data, geobox=dst_geobox, resampling=method) + reprojected = data.odc.reproject( + how=projection, resolution=resolution, resampling=method + ) - reprojected = data if "longitude" in reprojected.dims: reprojected = reprojected.rename({"longitude": "x"}) diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index 968845fc..cc169f27 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -28,74 +28,3 @@ def isnull(data): def notnull(data): return ~isnull(data) - - -class Unit(Enum): - metre = "metre" - degree = "degree" - - -def approx_metres_2_degrees(length_metres=10): - """Angle of rotation calculated off the earths equatorial radius in metres. - Ref: https://en.ans.wiki/5729/how-to-convert-meters-to-degrees/""" - earth_equator_radius_metres = 6378160 - return ((180 * length_metres) / (earth_equator_radius_metres * math.pi)) * 10 - - -def approx_degrees_2_metres(angle_in_degrees=0.0009): - """Calculated the metres travelled for an angle of rotation on the earths equatorial radius in metres. - Ref: https://en.ans.wiki/5728/how-to-convert-degrees-to-meters/""" - earth_equator_radius_metres = 6378160 - return ((earth_equator_radius_metres * math.pi * angle_in_degrees) / 180) / 10 - - -def detect_changing_unit(src_crs, dst_crs, src_res=None): - """Is there a unit change between the src and dst CRS. If so, return the approximate value for the new resolution.""" - - src_unit = Unit(src_crs.axis_info[0].unit_name) - dst_unit = Unit(dst_crs.axis_info[0].unit_name) - - # If not dst_res has been set, check the src_unit and dst_unit to - # determine whether the src_res needs to be converted. - if src_unit != dst_unit: - if dst_unit == Unit.metre: - print("deg 2 met") - return approx_degrees_2_metres(src_res) - if dst_unit == Unit.degree: - print("met 2 deg") - return approx_metres_2_degrees(src_res) - return src_res - - -def prepare_geobox( - data, dst_crs, dst_res, src_res, new_top_left_x, new_top_left_y, scale=False -): - """Get the destination geobox ready for the transformation.""" - # Docs for geotransform - # https://gdal.org/tutorials/geotransforms_tut.html - - new_affine = Affine( - dst_res, - 0, - new_top_left_x, - 0, - # Negative pixel width used for pixel height - -dst_res, - new_top_left_y, - ) - - if scale: - updated_resolution = resolution_from_affine(new_affine) - - x_min, x_max = min(data.x.values), max(data.x.values) - y_min, y_max = min(data.y.values), max(data.y.values) - - x_length = math.ceil(abs(x_min - x_max) / updated_resolution.x) - y_length = math.ceil(abs(y_min - y_max) / -updated_resolution.y) - - else: - x_length = len(data.x) - y_length = len(data.y) - - dst_geo = GeoBox(width=x_length, height=y_length, crs=dst_crs, affine=new_affine) - return dst_geo diff --git a/pyproject.toml b/pyproject.toml index 9a53dbd8..1c5c7d25 100644 --- a/pyproject.toml +++ b/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" From 830ee60aeca16f70fa0fd79a15c476a5b51c9e6f Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 25 May 2023 09:06:04 +0200 Subject: [PATCH 08/26] add reformat --- .../process_implementations/cubes/resample.py | 6 +++++- tests/test_resample.py | 1 + 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 3a28e9d3..e603e7af 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -41,7 +41,7 @@ def resample_spatial( f"[{', '.join(resample_methods_list)}]" ) - # Re-order, this is specifically done for xr_reproject + # Re-order, this is specifically done for odc reproject data = data.transpose("bands", "t", "y", "x") if not projection: @@ -66,4 +66,8 @@ def resample_spatial( reprojected = reprojected.rename({"latitude": "y"}) reprojected.attrs["crs"] = data.rio.crs + + # Undo odc specific re-ordering. + reprojected = reprojected.transpose("bands", "t", "x", "y") + return reprojected diff --git a/tests/test_resample.py b/tests/test_resample.py index 5ec4258d..35c9cdd1 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -44,6 +44,7 @@ def test_resample_spatial( verify_crs=False, ) + assert output_cube.odc.spatial_dims == ("y", "x") assert output_cube.rio.crs == CRS.from_user_input(output_crs) assert resolution_from_affine(output_cube.geobox.affine).x == output_res assert resolution_from_affine(output_cube.geobox.affine).y == -output_res From aeecb07f68e0b5695470712bf61073f4eee58d72 Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 25 May 2023 09:28:24 +0200 Subject: [PATCH 09/26] accept reformats --- .../process_implementations/cubes/resample.py | 6 ++++++ .../process_implementations/cubes/utils.py | 11 +++++++---- tests/test_resample.py | 7 +++++-- 3 files changed, 18 insertions(+), 6 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index e603e7af..461c07cd 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -4,6 +4,9 @@ from odc.geo.geobox import resolution_from_affine from pyproj.crs import CRS, CRSError +from openeo_processes_dask.process_implementations.cubes.utils import ( + approx_metres_2_degrees, +) from openeo_processes_dask.process_implementations.data_model import RasterCube resample_methods_list = [ @@ -55,6 +58,9 @@ def resample_spatial( if not resolution: resolution = resolution_from_affine(data.geobox.affine).x + if projection.is_geographic: + resolution = approx_metres_2_degrees(resolution) + reprojected = data.odc.reproject( how=projection, resolution=resolution, resampling=method ) diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index cc169f27..feb8733a 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -3,11 +3,7 @@ except ImportError: dask = None import math -from enum import Enum -from affine import Affine -from datacube.utils.geometry import GeoBox -from odc.geo.geobox import resolution_from_affine from xarray.core.duck_array_ops import isnull as xr_isnull @@ -28,3 +24,10 @@ def isnull(data): def notnull(data): return ~isnull(data) + + +def approx_metres_2_degrees(length_metres=10): + """Angle of rotation calculated off the earths equatorial radius in metres. + Ref: https://en.ans.wiki/5729/how-to-convert-meters-to-degrees/""" + earth_equator_radius_metres = 6378160 + return ((180 * length_metres) / (earth_equator_radius_metres * math.pi)) * 10 diff --git a/tests/test_resample.py b/tests/test_resample.py index 35c9cdd1..644652cf 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -16,6 +16,7 @@ 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]) @@ -46,5 +47,7 @@ def test_resample_spatial( assert output_cube.odc.spatial_dims == ("y", "x") assert output_cube.rio.crs == CRS.from_user_input(output_crs) - assert resolution_from_affine(output_cube.geobox.affine).x == output_res - assert resolution_from_affine(output_cube.geobox.affine).y == -output_res + + 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 From 5b034ed6c12613620fd0f2e67fe948144f8c0e8a Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 25 May 2023 09:41:45 +0200 Subject: [PATCH 10/26] accept reformatting --- .../process_implementations/cubes/resample.py | 3 --- .../process_implementations/cubes/utils.py | 7 ------- 2 files changed, 10 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 461c07cd..84d6efc8 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -58,9 +58,6 @@ def resample_spatial( if not resolution: resolution = resolution_from_affine(data.geobox.affine).x - if projection.is_geographic: - resolution = approx_metres_2_degrees(resolution) - reprojected = data.odc.reproject( how=projection, resolution=resolution, resampling=method ) diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index feb8733a..4ba8d994 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -24,10 +24,3 @@ def isnull(data): def notnull(data): return ~isnull(data) - - -def approx_metres_2_degrees(length_metres=10): - """Angle of rotation calculated off the earths equatorial radius in metres. - Ref: https://en.ans.wiki/5729/how-to-convert-meters-to-degrees/""" - earth_equator_radius_metres = 6378160 - return ((180 * length_metres) / (earth_equator_radius_metres * math.pi)) * 10 From 0fd1d14d2ebf827242bdf3027a26b5ad9c08c4fe Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 25 May 2023 09:44:42 +0200 Subject: [PATCH 11/26] drop broken import --- .../process_implementations/cubes/resample.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 84d6efc8..e603e7af 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -4,9 +4,6 @@ from odc.geo.geobox import resolution_from_affine from pyproj.crs import CRS, CRSError -from openeo_processes_dask.process_implementations.cubes.utils import ( - approx_metres_2_degrees, -) from openeo_processes_dask.process_implementations.data_model import RasterCube resample_methods_list = [ From 6b5496c11befd22f656f61b26fc17b7982b7d314 Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 25 May 2023 14:21:31 +0200 Subject: [PATCH 12/26] accept reformatting --- openeo_processes_dask/process_implementations/cubes/resample.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index e603e7af..5b6acdc5 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -29,7 +29,7 @@ def resample_spatial( method: str = "near", align: str = "upper-left", ): - """Resample the input rastercube by the provided epsg_code. Currently not accepting resolution resampling.""" + """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": From f5106544c8a1a0ac2e2e7bfe0c9b54ee19ae3474 Mon Sep 17 00:00:00 2001 From: sean Date: Mon, 5 Jun 2023 08:43:53 +0200 Subject: [PATCH 13/26] accept reformatting --- .../process_implementations/cubes/resample.py | 18 ++++++++++++------ tests/test_resample.py | 6 ++++++ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 5b6acdc5..c506b9ed 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -1,3 +1,4 @@ +import logging from typing import Union import odc.geo @@ -6,6 +7,8 @@ from openeo_processes_dask.process_implementations.data_model import RasterCube +logger = logging.getLogger(__name__) + resample_methods_list = [ "near", "bilinear", @@ -31,6 +34,9 @@ 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 align == "upper-left": + logging.warning("Warning: align parameter is unused by current implementation.") + # Assert resampling method is correct. if method == "near": method = "nearest" @@ -42,10 +48,10 @@ def resample_spatial( ) # Re-order, this is specifically done for odc reproject - data = data.transpose("bands", "t", "y", "x") + data_cp = data.transpose("bands", "t", "y", "x") if not projection: - projection = data.rio.crs + projection = data_cp.rio.crs try: projection = CRS(projection) @@ -53,9 +59,9 @@ def resample_spatial( raise CRSError(f"{projection} Can not be parsed to CRS.") if not resolution: - resolution = resolution_from_affine(data.geobox.affine).x + resolution = resolution_from_affine(data_cp.geobox.affine).x - reprojected = data.odc.reproject( + reprojected = data_cp.odc.reproject( how=projection, resolution=resolution, resampling=method ) @@ -65,9 +71,9 @@ def resample_spatial( if "latitude" in reprojected.dims: reprojected = reprojected.rename({"latitude": "y"}) - reprojected.attrs["crs"] = data.rio.crs + reprojected.attrs["crs"] = data_cp.rio.crs # Undo odc specific re-ordering. - reprojected = reprojected.transpose("bands", "t", "x", "y") + reprojected = reprojected.transpose(*data.dims) return reprojected diff --git a/tests/test_resample.py b/tests/test_resample.py index 644652cf..1554c7be 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -51,3 +51,9 @@ def test_resample_spatial( 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 + + assert min(output_cube.x) >= -180 + assert max(output_cube.x) <= 180 + + assert min(output_cube.y) >= -90 + assert max(output_cube.y) <= 90 From 4191013c08534ca47a77177cdd97275b130f96c5 Mon Sep 17 00:00:00 2001 From: sean Date: Mon, 5 Jun 2023 08:50:30 +0200 Subject: [PATCH 14/26] fixed broken tests --- tests/test_resample.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_resample.py b/tests/test_resample.py index 1554c7be..a1aa2edc 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -52,6 +52,7 @@ def test_resample_spatial( 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 From ea2efbeab1b6f663c0bcd732874cb6c1c9526c72 Mon Sep 17 00:00:00 2001 From: Lukas Weidenholzer Date: Tue, 6 Jun 2023 13:46:55 +0200 Subject: [PATCH 15/26] add comment --- .../process_implementations/cubes/resample.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index c506b9ed..332696ea 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -37,7 +37,8 @@ def resample_spatial( if align == "upper-left": logging.warning("Warning: align parameter is unused by current implementation.") - # Assert resampling method is correct. + # This is necessary, because the resampling methods in rasterio are used + # Reference here: https://rasterio.readthedocs.io/en/stable/api/rasterio.enums.html#rasterio.enums.Resampling if method == "near": method = "nearest" From f457719c211ef24d6cbd984e3625bfca8f31d6d3 Mon Sep 17 00:00:00 2001 From: Lukas Weidenholzer Date: Tue, 6 Jun 2023 13:47:08 +0200 Subject: [PATCH 16/26] add negative test case --- tests/test_resample.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/test_resample.py b/tests/test_resample.py index a1aa2edc..e2e9490c 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -34,6 +34,11 @@ def test_resample_spatial( 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 ) From a8e565be0cbefb96ed04aa8502466fcc5864e314 Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 22 Jun 2023 09:11:08 +0200 Subject: [PATCH 17/26] resolve locally --- .../process_implementations/cubes/resample.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index c506b9ed..62ad39d5 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -28,15 +28,11 @@ def resample_spatial( data: RasterCube, projection: Union[str, int] = None, - resolution: Union[str, int] = None, + resolution: int = 0, 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" @@ -48,7 +44,12 @@ def resample_spatial( ) # Re-order, this is specifically done for odc reproject - data_cp = data.transpose("bands", "t", "y", "x") + data_cp = data.transpose( + data.openeo.band_dims[0], + data.openeo.temporal_dims[0], + data.openeo.y_dim, + data.openeo.x_dim, + ) if not projection: projection = data_cp.rio.crs From 4473660d2fabcf4d30551fc497dd0b0a94349fe3 Mon Sep 17 00:00:00 2001 From: sean Date: Fri, 30 Jun 2023 10:01:10 +0200 Subject: [PATCH 18/26] try merge from main to fix tests and change crs function --- openeo_processes_dask/process_implementations/cubes/resample.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 62ad39d5..915ef909 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -55,7 +55,7 @@ def resample_spatial( projection = data_cp.rio.crs try: - projection = CRS(projection) + projection = CRS.from_user_input(projection) except CRSError as e: raise CRSError(f"{projection} Can not be parsed to CRS.") From c3c06383508fa0fd1aacbef38d37a04c1d15018e Mon Sep 17 00:00:00 2001 From: Sean <63349506+SerRichard@users.noreply.github.com> Date: Thu, 17 Aug 2023 09:02:08 +0200 Subject: [PATCH 19/26] Update openeo_processes_dask/process_implementations/cubes/resample.py Co-authored-by: Lukas Weidenholzer <17790923+LukeWeidenwalker@users.noreply.github.com> --- openeo_processes_dask/process_implementations/cubes/resample.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 6751bea5..9d120545 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -27,7 +27,7 @@ def resample_spatial( data: RasterCube, - projection: Union[str, int] = None, + projection: Optional[Union[str, int]] = None, resolution: int = 0, method: str = "near", ): From dc9a4ff794f97c648b306a2aefb05bb144ca7fb0 Mon Sep 17 00:00:00 2001 From: Sean <63349506+SerRichard@users.noreply.github.com> Date: Thu, 17 Aug 2023 09:02:19 +0200 Subject: [PATCH 20/26] Update openeo_processes_dask/process_implementations/cubes/resample.py Co-authored-by: Lukas Weidenholzer <17790923+LukeWeidenwalker@users.noreply.github.com> --- openeo_processes_dask/process_implementations/cubes/resample.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 9d120545..47901980 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -51,7 +51,7 @@ def resample_spatial( data.openeo.x_dim, ) - if not projection: + if projection is None: projection = data_cp.rio.crs try: From 6ce40e738e53c252565952d0c8c0b1cff4bdb17e Mon Sep 17 00:00:00 2001 From: Sean <63349506+SerRichard@users.noreply.github.com> Date: Thu, 17 Aug 2023 09:02:27 +0200 Subject: [PATCH 21/26] Update openeo_processes_dask/process_implementations/cubes/resample.py Co-authored-by: Lukas Weidenholzer <17790923+LukeWeidenwalker@users.noreply.github.com> --- openeo_processes_dask/process_implementations/cubes/resample.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 47901980..3a1cdf10 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -59,7 +59,7 @@ def resample_spatial( except CRSError as e: raise CRSError(f"{projection} Can not be parsed to CRS.") - if not resolution: + if resolution == 0: resolution = resolution_from_affine(data_cp.odc.geobox.affine).x reprojected = data_cp.odc.reproject( From 4132feafb63f60294c5902c75696d11a87959de8 Mon Sep 17 00:00:00 2001 From: Sean <63349506+SerRichard@users.noreply.github.com> Date: Thu, 17 Aug 2023 09:02:52 +0200 Subject: [PATCH 22/26] Update openeo_processes_dask/process_implementations/cubes/resample.py Co-authored-by: Lukas Weidenholzer <17790923+LukeWeidenwalker@users.noreply.github.com> --- openeo_processes_dask/process_implementations/cubes/resample.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 3a1cdf10..798ff52f 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -57,7 +57,7 @@ def resample_spatial( try: projection = CRS.from_user_input(projection) except CRSError as e: - raise CRSError(f"{projection} Can not be parsed to CRS.") + 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 From 251f8e8b757b01656ce5e5d4fc7be0dc4a397dfe Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 17 Aug 2023 09:05:56 +0200 Subject: [PATCH 23/26] comments updates --- .../process_implementations/cubes/resample.py | 7 ++++--- .../process_implementations/cubes/utils.py | 1 - 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 6751bea5..62d446a4 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -6,6 +6,7 @@ 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__) @@ -38,7 +39,7 @@ def resample_spatial( method = "nearest" elif method not in resample_methods_list: - raise Exception( + raise OpenEOException( f'Selected resampling method "{method}" is not available! Please select one of ' f"[{', '.join(resample_methods_list)}]" ) @@ -66,10 +67,10 @@ def resample_spatial( how=projection, resolution=resolution, resampling=method ) - if "longitude" in reprojected.dims: + if "longitude" in reprojected.dims and "longitude" in data.dims: reprojected = reprojected.rename({"longitude": "x"}) - if "latitude" in reprojected.dims: + if "latitude" in reprojected.dims and "latitude" in data.dims: reprojected = reprojected.rename({"latitude": "y"}) reprojected.attrs["crs"] = data_cp.rio.crs diff --git a/openeo_processes_dask/process_implementations/cubes/utils.py b/openeo_processes_dask/process_implementations/cubes/utils.py index 4ba8d994..1bccf1bb 100644 --- a/openeo_processes_dask/process_implementations/cubes/utils.py +++ b/openeo_processes_dask/process_implementations/cubes/utils.py @@ -2,7 +2,6 @@ import dask except ImportError: dask = None -import math from xarray.core.duck_array_ops import isnull as xr_isnull From 3e08f66c54193e9aae34cc2b9f029fb936167e51 Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 17 Aug 2023 09:07:10 +0200 Subject: [PATCH 24/26] import typing optional --- .../process_implementations/cubes/resample.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index d92dca2f..5b9eced0 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -1,5 +1,5 @@ import logging -from typing import Union +from typing import Optional, Union import odc.geo.xr from odc.geo.geobox import resolution_from_affine @@ -58,7 +58,9 @@ def resample_spatial( 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 + 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 From 8dff319952ebf5e8c38eb7528dcc238952621598 Mon Sep 17 00:00:00 2001 From: sean Date: Thu, 17 Aug 2023 09:23:11 +0200 Subject: [PATCH 25/26] fix tessts --- .../process_implementations/cubes/resample.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py index 5b9eced0..01d4bb29 100644 --- a/openeo_processes_dask/process_implementations/cubes/resample.py +++ b/openeo_processes_dask/process_implementations/cubes/resample.py @@ -69,10 +69,10 @@ def resample_spatial( how=projection, resolution=resolution, resampling=method ) - if "longitude" in reprojected.dims and "longitude" in data.dims: + if "longitude" in reprojected.dims and "x" in data.dims: reprojected = reprojected.rename({"longitude": "x"}) - if "latitude" in reprojected.dims and "latitude" in data.dims: + if "latitude" in reprojected.dims and "y" in data.dims: reprojected = reprojected.rename({"latitude": "y"}) reprojected.attrs["crs"] = data_cp.rio.crs From 57af9bee6f07a114d7060381463d205c76b67575 Mon Sep 17 00:00:00 2001 From: Lukas Weidenholzer Date: Thu, 17 Aug 2023 11:11:13 +0200 Subject: [PATCH 26/26] bump process specs --- openeo_processes_dask/specs/openeo-processes | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openeo_processes_dask/specs/openeo-processes b/openeo_processes_dask/specs/openeo-processes index 730c13c5..51cdcd8b 160000 --- a/openeo_processes_dask/specs/openeo-processes +++ b/openeo_processes_dask/specs/openeo-processes @@ -1 +1 @@ -Subproject commit 730c13c50fa9b0e3bd9525f950dc82bbe2799329 +Subproject commit 51cdcd8b8bb79cf5caf89222b9c2ac8a4effb89b