Add utils.transform_resolution

CHANGES.md

@@ -1,3 +1,10 @@
+## Changes in 0.3.4
+
+- Added `utils.transform_resolution` to convert spatial resolution between coordinate 
+  reference systems (CRS) and introduced deprecation warnings for 
+  `utils.resolution_degrees_to_meters` and `utils.resolution_meters_to_degrees`, which
+  are superseded by `utils.transform_resolution`.
+
 ## Changes in 0.3.3
 
 - Added function `utils.resolution_degrees_to_meters` to convert spatial resolution

tests/test_utils.py

@@ -30,6 +30,7 @@
     reproject_bbox,
     resolution_degrees_to_meters,
     resolution_meters_to_degrees,
+    transform_resolution,
 )
 
 from .sampledata import create_2x4x4_dataset_with_irregular_coords
@@ -418,9 +419,34 @@ def test_bbox_overlap(self):
         target = (-170, 0, -160, 10)
         self.assertAlmostEqual(0.0, bbox_overlap(source, target))
 
+    def test_transform_resolution(self):
+        # EPSG:4326 → UTM (meters)
+        res = transform_resolution((3, 60), 1.0, "EPSG:4326", "EPSG:32631")
+        self.assertAlmostEqual(55777.9, res[0], places=1)
+        self.assertAlmostEqual(111376.2, res[1], places=1)
+
+        # UTM (meters) → EPSG:4326
+        res = transform_resolution(
+            (500_000, 6_651_411), 111376, "EPSG:32631", "EPSG:4326"
+        )
+        self.assertAlmostEqual(2.0, res[0], places=2)
+        self.assertAlmostEqual(1.0, res[1], places=2)
+
+        # UTM → UTM (same CRS, should remain ~unchanged in magnitude)
+        res = transform_resolution(
+            (500000, 0), (1000, 1000), "EPSG:32631", "EPSG:32631"
+        )
+        self.assertAlmostEqual(1000, res[0])
+        self.assertAlmostEqual(1000, res[1])
+
+        # US survey foot to meter
+        res = transform_resolution((0, 0), (1, 1), "EPSG:3561", "EPSG:32605")
+        self.assertAlmostEqual(0.30525, res[0], places=4)
+        self.assertAlmostEqual(0.30525, res[1], places=4)
+
     def test_resolution_meters_to_degrees(self):
         # 111320 m ≈ 1 degree at equator
-        lat_deg, lon_deg = resolution_meters_to_degrees(111320, 0)
+        lon_deg, lat_deg = resolution_meters_to_degrees(111320, 0)
         self.assertAlmostEqual(1.0, lat_deg, places=6)
         self.assertAlmostEqual(1.0, lon_deg, places=6)
 
@@ -436,7 +462,7 @@ def test_resolution_meters_to_degrees(self):
 
     def test_resolution_degrees_to_meters(self):
         # 111320 m ≈ 1 degree at equator
-        lat_deg, lon_deg = resolution_degrees_to_meters(1, 0)
+        lon_deg, lat_deg = resolution_degrees_to_meters(1, 0)
         self.assertAlmostEqual(111320, lat_deg, places=6)
         self.assertAlmostEqual(111320, lon_deg, places=6)
 

xcube_resampling/affine.py

@@ -22,7 +22,6 @@
 import math
 from collections.abc import Iterable, Sequence
 
-import dask
 import dask.array as da
 import numpy as np
 import xarray as xr

xcube_resampling/utils.py

@@ -19,6 +19,7 @@
 # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 # DEALINGS IN THE SOFTWARE.
 
+import warnings
 from collections.abc import Callable, Hashable, Iterable, Mapping, Sequence
 from dataclasses import dataclass
 
@@ -313,12 +314,54 @@ def _split_bbox_antimeridian(bbox: Sequence[FloatInt]) -> Sequence[Sequence[Floa
     return [bbox]
 
 
+def transform_resolution(
+    ref_point: tuple[FloatInt, FloatInt],
+    resolution: FloatInt | tuple[FloatInt, FloatInt],
+    src_crs: str | pyproj.CRS,
+    dst_crs: str | pyproj.CRS,
+) -> tuple[FloatInt, FloatInt]:
+    """Estimate local spatial resolution in the destination CRS using
+    finite differences.
+
+    Args:
+        ref_point: Reference point (easting, northing) in the source CRS.
+        resolution: Spatial resolution in the source CRS. Can be a single number
+            (applied equally to both axes) or a tuple `(x_res, y_res)`.
+        src_crs: Source coordinate reference system.
+        dst_crs: Destination coordinate reference system.
+
+    Returns:
+        tuple: Estimated local resolution in the destination CRS as `(x_res, y_res)`.
+    """
+    transformer = pyproj.Transformer.from_crs(src_crs, dst_crs, always_xy=True)
+    if not isinstance(resolution, tuple):
+        resolution = (resolution, resolution)
+
+    # reference point
+    x0, y0 = transformer.transform(*ref_point)
+
+    # step in x direction
+    x1, y1 = transformer.transform(ref_point[0] + resolution[0], ref_point[1])
+
+    # step in y direction
+    x2, y2 = transformer.transform(ref_point[0], ref_point[1] + resolution[1])
+
+    # Euclidean distances
+    res_x = np.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2)
+    res_y = np.sqrt((x2 - x0) ** 2 + (y2 - y0) ** 2)
+
+    return res_x, res_y
+
+
 def resolution_meters_to_degrees(
     resolution: FloatInt | tuple[FloatInt, FloatInt], latitude: FloatInt
 ) -> tuple[FloatInt, FloatInt]:
     """Convert spatial resolution from meters to degrees in longitude and latitude
     at a given geographic latitude.
 
+    This function is deprecated and will be removed in a future release.
+    Use `transform_resolution` instead for CRS-aware and more accurate conversions.
+
     Args:
         resolution: Spatial resolution in meters. Can be a single number
             (applied equally to both axes) or a tuple ``(x_res, y_res)``.
@@ -333,6 +376,13 @@ def resolution_meters_to_degrees(
         - 1 degree of longitude ≈ 111,320 * cos(latitude) meters.
 
     """
+    warnings.warn(
+        "resolution_meters_to_degrees is deprecated and will be removed in a future release. "
+        "Use `transform_resolution` instead for accurate CRS-based resolution handling.",
+        DeprecationWarning,
+        stacklevel=2,
+    )
+
     if not isinstance(resolution, tuple):
         resolution = (resolution, resolution)
     return (
@@ -346,6 +396,9 @@ def resolution_degrees_to_meters(
 ) -> tuple[FloatInt, FloatInt]:
     """Convert spatial resolution from degrees to meters at a given geographic latitude.
 
+    This function is deprecated and will be removed in a future release.
+    Use `transform_resolution` instead for CRS-aware and more accurate conversions.
+
     Args:
         resolution: Spatial resolution in degrees. Can be a single number
             (applied equally to both axes) or a tuple ``(lon_res, lat_res)``.
@@ -360,6 +413,13 @@ def resolution_degrees_to_meters(
         - 1 degree of longitude ≈ 111,320 * cos(latitude) meters.
 
     """
+    warnings.warn(
+        "resolution_degrees_to_meters is deprecated and will be removed in a future release. "
+        "Use `transform_resolution` instead for accurate CRS-based resolution handling.",
+        DeprecationWarning,
+        stacklevel=2,
+    )
+
     if not isinstance(resolution, tuple):
         resolution = (resolution, resolution)
     return (

xcube_resampling/version.py

@@ -19,4 +19,4 @@
 # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 # DEALINGS IN THE SOFTWARE.
 
-__version__ = "0.3.3"
+__version__ = "0.3.4"