From 7bbed7493a8d23b5ec33a6888065707e363660f6 Mon Sep 17 00:00:00 2001 From: "Gregory H. Halverson" Date: Tue, 19 Aug 2025 12:52:21 -0700 Subject: [PATCH 1/2] fixing transformations --- rasters/raster_geometry.py | 13 ++------- rasters/transform_xy.py | 54 ++++++++++++++++++++++++++++++++++++ tests/test_transform_xy.py | 56 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 113 insertions(+), 10 deletions(-) create mode 100644 rasters/transform_xy.py create mode 100644 tests/test_transform_xy.py diff --git a/rasters/raster_geometry.py b/rasters/raster_geometry.py index bd0e48b..13fe45a 100644 --- a/rasters/raster_geometry.py +++ b/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 diff --git a/rasters/transform_xy.py b/rasters/transform_xy.py new file mode 100644 index 0000000..5d7ad0d --- /dev/null +++ b/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 diff --git a/tests/test_transform_xy.py b/tests/test_transform_xy.py new file mode 100644 index 0000000..d4085fc --- /dev/null +++ b/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 From 64dc94af623c423e3e97ec52fdfa714d24ba930f Mon Sep 17 00:00:00 2001 From: "Gregory H. Halverson" Date: Tue, 19 Aug 2025 12:52:44 -0700 Subject: [PATCH 2/2] v1.11.0 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index cdad46a..d065729 100644 --- a/pyproject.toml +++ b/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 = [