Merge pull request #51 from gregory-halverson-jpl/main

fixing transformations

Author: gregory-halverson-jpl

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "rasters"
-version = "1.10.0"
+version = "1.11.0"
 description = "raster processing toolkit"
 readme = "README.md"
 authors = [

rasters/raster_geometry.py

@@ -590,9 +590,10 @@ def geolocation(self) -> RasterGeolocation:
 
         return RasterGeolocation(x=self.x, y=self.y, crs=self.crs)
 
-    def to_crs(self, crs: Union[CRS, str] = WGS84) -> RasterGeolocation:
+    def to_crs(self, crs: Union[CRS, str] = WGS84) -> 'RasterGeolocation':
         from .CRS import CRS
         from .raster_geolocation import RasterGeolocation
+        from .transform_xy import transform_xy
 
         # validate destination CRS
         if not isinstance(crs, CRS):
@@ -601,15 +602,7 @@ def to_crs(self, crs: Union[CRS, str] = WGS84) -> RasterGeolocation:
         if self.crs == crs:
             return self
 
-        if self.is_geographic:
-            x, y = shapely_transform(self.crs, crs, self.y, self.x)
-        else:
-            x, y = shapely_transform(self.crs, crs, self.x, self.y)
-
-        if crs.is_geographic:
-            x = np.where(np.logical_or(x < -180, x > 180), np.nan, x)
-            y = np.where(np.logical_or(y < -90, y > 90), np.nan, y)
-
+        x, y = transform_xy(x=self.x, y=self.y, source_crs=self.crs, target_crs=crs)
         geolocation = RasterGeolocation(x=x, y=y, crs=crs)
 
         return geolocation

rasters/transform_xy.py

@@ -0,0 +1,54 @@
+from pyproj import Transformer
+import numpy as np
+
+def transform_xy(
+    x: np.ndarray,
+    y: np.ndarray,
+    source_crs,
+    target_crs
+) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Transform arrays of x and y coordinates from a source CRS to a target CRS.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        Array of x coordinates (e.g., longitude or easting).
+    y : np.ndarray
+        Array of y coordinates (e.g., latitude or northing).
+    source_crs : pyproj.CRS, str, or custom CRS
+        The source coordinate reference system. Can be a pyproj.CRS, proj4 string, EPSG code, or custom CRS class.
+    target_crs : pyproj.CRS, str, or custom CRS
+        The target coordinate reference system. Can be a pyproj.CRS, proj4 string, EPSG code, or custom CRS class.
+
+    Returns
+    -------
+    x_t : np.ndarray
+        Transformed x coordinates in the target CRS.
+    y_t : np.ndarray
+        Transformed y coordinates in the target CRS.
+
+    Notes
+    -----
+    - Handles both geographic and projected CRS.
+    - If the target CRS is geographic, output coordinates outside valid bounds are set to np.nan.
+    - Accepts pyproj.CRS, string, or custom CRS class with 'is_geographic' attribute.
+    """
+    try:
+        from .CRS import CRS
+        if not isinstance(source_crs, CRS):
+            source_crs = CRS(source_crs)
+        if not isinstance(target_crs, CRS):
+            target_crs = CRS(target_crs)
+    except ImportError:
+        pass
+
+    transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)
+    x_t, y_t = transformer.transform(x, y)
+
+    # If target CRS is geographic, clip to valid bounds
+    if hasattr(target_crs, 'is_geographic') and target_crs.is_geographic:
+        x_t = np.where(np.logical_or(x_t < -180, x_t > 180), np.nan, x_t)
+        y_t = np.where(np.logical_or(y_t < -90, y_t > 90), np.nan, y_t)
+
+    return x_t, y_t

tests/test_transform_xy.py

@@ -0,0 +1,56 @@
+import numpy as np
+from pyproj import CRS
+from rasters.transform_xy import transform_xy
+
+def test_transform_xy_geographic_to_projected():
+    # WGS84 geographic CRS
+    source_crs = CRS.from_epsg(4326)
+    # UTM zone 33N projected CRS
+    target_crs = CRS.from_epsg(32633)
+    # Example coordinates: longitude, latitude
+    x = np.array([12.0, 13.0])
+    y = np.array([55.0, 56.0])
+    x_t, y_t = transform_xy(x, y, source_crs, target_crs)
+    # Check output shape and type
+    assert x_t.shape == x.shape
+    assert y_t.shape == y.shape
+    assert np.all(np.isfinite(x_t))
+    assert np.all(np.isfinite(y_t))
+    # Check that transformed coordinates are not equal to input
+    assert not np.allclose(x, x_t)
+    assert not np.allclose(y, y_t)
+
+def test_transform_xy_projected_to_geographic():
+    # UTM zone 33N projected CRS
+    source_crs = CRS.from_epsg(32633)
+    # WGS84 geographic CRS
+    target_crs = CRS.from_epsg(4326)
+    # Example coordinates: easting, northing
+    x = np.array([500000, 600000])
+    y = np.array([6100000, 6200000])
+    x_t, y_t = transform_xy(x, y, source_crs, target_crs)
+    # Check output shape and type
+    assert x_t.shape == x.shape
+    assert y_t.shape == y.shape
+    assert np.all(np.isfinite(x_t))
+    assert np.all(np.isfinite(y_t))
+    # Check that transformed coordinates are not equal to input
+    assert not np.allclose(x, x_t)
+    assert not np.allclose(y, y_t)
+
+def test_transform_xy_clip_geographic():
+    # WGS84 geographic CRS
+    source_crs = CRS.from_epsg(4326)
+    target_crs = CRS.from_epsg(4326)
+    # Coordinates outside valid bounds
+    x = np.array([-200, 0, 200])
+    y = np.array([-100, 0, 100])
+    x_t, y_t = transform_xy(x, y, source_crs, target_crs)
+    # Out-of-bounds should be nan
+    assert np.isnan(x_t[0])
+    assert np.isnan(x_t[2])
+    assert np.isnan(y_t[0])
+    assert np.isnan(y_t[2])
+    # Valid should remain
+    assert x_t[1] == 0
+    assert y_t[1] == 0