From e1c91466c4c477159ae548d1a06ffe4d661eb27a Mon Sep 17 00:00:00 2001 From: "Gregory H. Halverson" Date: Tue, 19 Aug 2025 10:45:03 -0700 Subject: [PATCH 1/2] unit test for writing GeoTIFF --- tests/test_write_geotiff.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/test_write_geotiff.py diff --git a/tests/test_write_geotiff.py b/tests/test_write_geotiff.py new file mode 100644 index 0000000..e69de29 From 01dae6b9544283a10d5ea0844fbeb8e33dad5b4c Mon Sep 17 00:00:00 2001 From: "Gregory H. Halverson" Date: Tue, 19 Aug 2025 11:26:36 -0700 Subject: [PATCH 2/2] fixing UTM projection --- pyproject.toml | 2 +- rasters/local_UTM_proj4.py | 22 ++-- rasters/local_UTM_proj4_from_lat_lon.py | 19 ---- rasters/raster_geometry.py | 17 +-- rasters/spatial_geometry.py | 27 +---- tests/test_local_UTM_proj4.py | 46 ++++++++ tests/test_raster_geometry_local_UTM_proj4.py | 102 ++++++++++++++++++ 7 files changed, 176 insertions(+), 59 deletions(-) create mode 100644 tests/test_local_UTM_proj4.py create mode 100644 tests/test_raster_geometry_local_UTM_proj4.py diff --git a/pyproject.toml b/pyproject.toml index 0760f7f..cdad46a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "rasters" -version = "1.9.0" +version = "1.10.0" description = "raster processing toolkit" readme = "README.md" authors = [ diff --git a/rasters/local_UTM_proj4.py b/rasters/local_UTM_proj4.py index 86bf21f..23891e2 100644 --- a/rasters/local_UTM_proj4.py +++ b/rasters/local_UTM_proj4.py @@ -6,12 +6,12 @@ from .CRS import CRS from .point import Point -def local_UTM_proj4(cls, point_latlon: Union[Point, str]) -> CRS: +def local_UTM_proj4(point_latlon: Union[Point, str]) -> CRS: """ Generate a local UTM projection based on a given point. Args: - point_latlon (Union[Point, str]): Point object or WKT string containing latitude and longitude. + point_latlon (Union[Point, str]): Point object or WKT string containing latitude and longitude. Returns: CRS: pyproj.CRS object of the local UTM projection. @@ -28,9 +28,17 @@ def local_UTM_proj4(cls, point_latlon: Union[Point, str]) -> CRS: if not (-180 <= lon <= 180): raise ValueError("Longitude must be between -180 and 180 degrees") - UTM_zone = (math.floor((lon + 180) / 6) % 60) + 1 - UTM_proj4 = CRS( - f"+proj=utm +zone={UTM_zone} {'+south ' if lat < 0 else ''}+ellps=WGS84 +datum=WGS84 +units=m +no_defs" - ) - + UTM_zone = int(math.floor((lon + 180) / 6)) + UTM_zone = max(1, min(60, UTM_zone)) # Clamp to valid UTM zones + import warnings + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="You will likely lose important projection information when converting to a PROJ string from another format.", + category=UserWarning, + module="pyproj" + ) + UTM_proj4 = CRS( + f"+proj=utm +zone={UTM_zone} {'+south ' if lat < 0 else ''}+ellps=WGS84 +datum=WGS84 +units=m +no_defs" + ) return UTM_proj4 diff --git a/rasters/local_UTM_proj4_from_lat_lon.py b/rasters/local_UTM_proj4_from_lat_lon.py index 2b9de4b..10de6a0 100644 --- a/rasters/local_UTM_proj4_from_lat_lon.py +++ b/rasters/local_UTM_proj4_from_lat_lon.py @@ -1,20 +1 @@ -import numpy as np - def local_UTM_proj4_from_lat_lon(lat: float, lon: float) -> str: - """ - Generate a local UTM projection string (proj4 format) based on given latitude and longitude. - - This function calculates the UTM zone based on the provided longitude and - constructs a proj4 string for the corresponding UTM projection. - - Args: - lat (float): The latitude of the point. - lon (float): The longitude of the point. - - Returns: - str: The proj4 string for the local UTM projection. - """ - UTM_zone = (np.floor((lon + 180) / 6) % 60) + 1 - UTM_proj4 = f"+proj=utm +zone={int(UTM_zone)} {'+south ' if lat < 0 else ''}+ellps=WGS84 +datum=WGS84 +units=m +no_defs" - - return UTM_proj4 diff --git a/rasters/raster_geometry.py b/rasters/raster_geometry.py index 68c7d7d..bd0e48b 100644 --- a/rasters/raster_geometry.py +++ b/rasters/raster_geometry.py @@ -13,6 +13,7 @@ from scipy.ndimage import shift from shapely.geometry.base import BaseGeometry from shapely.ops import transform as shapely_transform +from .local_UTM_proj4 import local_UTM_proj4 from .CRS import WGS84 from .spatial_geometry import SpatialGeometry @@ -231,13 +232,15 @@ def local_UTM_EPSG(self) -> str: @property def local_UTM_proj4(self) -> str: - centroid = self.centroid.latlon - lat = centroid.y - lon = centroid.x - UTM_zone = (np.floor((lon + 180) / 6) % 60) + 1 - UTM_proj4 = f"+proj=utm +zone={UTM_zone} {'+south ' if lat < 0 else ''}+ellps=WGS84 +datum=WGS84 +units=m +no_defs" - - return UTM_proj4 + import warnings + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="You will likely lose important projection information when converting to a PROJ string from another format.", + category=UserWarning, + module="pyproj" + ) + return local_UTM_proj4(self.centroid.latlon).to_proj4() @property def is_point(self) -> bool: diff --git a/rasters/spatial_geometry.py b/rasters/spatial_geometry.py index fbf1e8a..70a5926 100644 --- a/rasters/spatial_geometry.py +++ b/rasters/spatial_geometry.py @@ -4,7 +4,6 @@ from typing import Union, TYPE_CHECKING from .CRS import CRS, WGS84 -from .local_UTM_proj4_from_lat_lon import local_UTM_proj4_from_lat_lon if TYPE_CHECKING: from .bbox import BBox @@ -119,28 +118,6 @@ def centroid_latlon(self) -> "Point": @property def local_UTM_proj4(self) -> str: - """ - Gets the proj4 string for the local UTM zone of the geometry. - - The local UTM zone is determined by the centroid of the geometry in the - WGS84 coordinate reference system (CRS). - - Returns: - str: The proj4 string for the local UTM zone of the geometry. - """ - + from .local_UTM_proj4 import local_UTM_proj4 centroid = self.centroid.latlon - lat = centroid.y - lon = centroid.x - - # Calculate the UTM zone using the centroid's longitude. - # UTM_zone = (np.floor((lon + 180) / 6) % 60) + 1 - - # Construct the proj4 string using the calculated UTM zone and the - # centroid's latitude to determine the hemisphere. - # UTM_proj4 = f"+proj=utm +zone={UTM_zone} {'+south ' if lat < 0 else ''}+ellps=WGS84 +datum=WGS84 +units=m +no_defs" - - # Use the helper function to get the local UTM proj4 string - UTM_proj4 = local_UTM_proj4_from_lat_lon(lat, lon) - - return UTM_proj4 + return local_UTM_proj4(centroid) diff --git a/tests/test_local_UTM_proj4.py b/tests/test_local_UTM_proj4.py new file mode 100644 index 0000000..4405d2e --- /dev/null +++ b/tests/test_local_UTM_proj4.py @@ -0,0 +1,46 @@ +import unittest +import warnings +from rasters.local_UTM_proj4 import local_UTM_proj4 +from rasters.point import Point +from pyproj import CRS + +class TestLocalUTMProj4(unittest.TestCase): + def test_utm_zone_integer(self): + # Test for a point in UTM zone 10 + point = Point(-120.0, 40.0) # lon, lat + crs = local_UTM_proj4(point) + self.assertIsInstance(crs, CRS) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="You will likely lose important projection information when converting to a PROJ string from another format.", + category=UserWarning, + module="pyproj" + ) + proj_str = crs.to_proj4() + self.assertIn("+zone=10 ", proj_str) + self.assertIn("+proj=utm", proj_str) + self.assertIn("+datum=WGS84", proj_str) + + def test_utm_zone_southern_hemisphere(self): + # Test for a point in southern hemisphere + point = Point(30.0, -20.0) # lon, lat + crs = local_UTM_proj4(point) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="You will likely lose important projection information when converting to a PROJ string from another format.", + category=UserWarning, + module="pyproj" + ) + proj_str = crs.to_proj4() + self.assertIn("+south", proj_str) + + def test_invalid_longitude(self): + # Test for invalid longitude + point = Point(200.0, 40.0) + with self.assertRaises(ValueError): + local_UTM_proj4(point) + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_raster_geometry_local_UTM_proj4.py b/tests/test_raster_geometry_local_UTM_proj4.py new file mode 100644 index 0000000..44da06e --- /dev/null +++ b/tests/test_raster_geometry_local_UTM_proj4.py @@ -0,0 +1,102 @@ + +import unittest +from rasters.raster_geometry import RasterGeometry +from rasters.point import Point +from rasters.CRS import CRS, WGS84 + +class DummyRasterGeometry(RasterGeometry): + def __init__(self, x_min, x_max, y_min, y_max, crs=WGS84): + super().__init__(crs=crs) + self._x_min = x_min + self._x_max = x_max + self._y_min = y_min + self._y_max = y_max + self._rows = 1 + self._cols = 1 + + @property + def x(self): + return [[self._x_min]] + + @property + def y(self): + return [[self._y_min]] + + @property + def xy(self): + return (self.x, self.y) + + @property + def rows(self): + return self._rows + + @property + def cols(self): + return self._cols + + @property + def cell_width(self): + return self._x_max - self._x_min + + @property + def cell_height(self): + return self._y_max - self._y_min + + @property + def x_min(self): + return self._x_min + + @property + def x_max(self): + return self._x_max + + @property + def y_min(self): + return self._y_min + + @property + def y_max(self): + return self._y_max + + @property + def grid(self): + return None + + @property + def corner_polygon(self): + return None + + @property + def boundary(self): + return None + + def _subset_index(self, y_slice, x_slice): + return self + + def _slice_coords(self, y_slice, x_slice): + return self + + def index_point(self, point): + return (0, 0) + + def index(self, geometry): + return (0, 0) + + def resize(self, dimensions): + return self + + def to_dict(self): + return {} + +class TestLocalUTMProj4Property(unittest.TestCase): + def test_local_UTM_proj4_property(self): + # Create a dummy raster centered at lon=-120, lat=40 + raster = DummyRasterGeometry(x_min=-120.0, x_max=-120.0, y_min=40.0, y_max=40.0) + proj4_str = raster.local_UTM_proj4 + self.assertIsInstance(proj4_str, str) + self.assertIn("+proj=utm", proj4_str) + self.assertIn("+zone=10", proj4_str) + self.assertIn("+datum=WGS84", proj4_str) + +if __name__ == "__main__": + unittest.main()