fixing UTM projection

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 = [

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

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

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:

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)

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()

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()