Replace remaining TransformTo (#101)

CHANGES.rst

@@ -4,7 +4,10 @@ Changelog of dask-geomodeling
 2.3.12 (unreleased)
 -------------------
 
-- Nothing changed yet.
+- Perform more geometry transformations via pyproj (see 2.3.11 notes).
+
+- The .geometry property of elementwise, reduction and combine RasterBlocks is
+  now computed through the extents of the arguments (and will always be a box). 
 
 
 2.3.11 (2022-11-22)

dask_geomodeling/geometry/aggregate.py

@@ -281,7 +281,7 @@ def get_sources_and_requests(self, **request):
             ]
 
         # transform the extent into the projection in which we aggregate
-        x1, y1, x2, y2 = utils.transform_extent(extent, req_srs, agg_srs)
+        x1, y1, x2, y2 = utils.Extent(extent, req_srs).transformed(agg_srs).bbox
 
         # estimate the amount of required pixels
         required_pixels = int(((x2 - x1) * (y2 - y1)) / (self.pixel_size ** 2))

dask_geomodeling/geometry/constructive.py

@@ -3,7 +3,7 @@
 """
 import numbers
 
-from dask_geomodeling.utils import Extent, shapely_transform, transform_extent
+from dask_geomodeling.utils import Extent, shapely_transform
 
 from .base import BaseSingle
 
@@ -83,10 +83,8 @@ def process(data, kwargs):
             req_srs = data["projection"]
             buf_srs = kwargs["buf_srs"]
             distance = kwargs["distance"]
-            extent = transform_extent(data["extent"], req_srs, buf_srs)
-            extent = Extent(extent, buf_srs).buffered(distance).bbox
-            extent = transform_extent(extent, buf_srs, req_srs)
-            return {"extent": extent, "projection": req_srs}
+            extent = Extent(data["extent"], req_srs).transformed(buf_srs).buffered(distance).transformed(req_srs)
+            return {"extent": extent.bbox, "projection": req_srs}
         else:
             raise NotImplementedError("Dunno this mode!")
 

dask_geomodeling/raster/combine.py

@@ -6,7 +6,7 @@
 import numpy as np
 
 from dask_geomodeling.utils import filter_none, get_dtype_max, get_index
-from dask_geomodeling.utils import GeoTransform
+from dask_geomodeling.utils import GeoTransform, Extent
 
 from .base import RasterBlock
 
@@ -107,14 +107,10 @@ def geometry(self):
             return
         elif len(geometries) == 1:
             return geometries[0]
-        result = geometries[0]
-        sr = result.GetSpatialReference()
+        extent = Extent.from_geometry(geometries[0])
         for geometry in geometries[1:]:
-            if not geometry.GetSpatialReference().IsSame(sr):
-                geometry = geometry.Clone()
-                geometry.TransformTo(sr)
-            result = result.Union(geometry)
-        return result
+            extent = extent.union(Extent.from_geometry(geometry))
+        return extent.as_geometry()
 
     @property
     def projection(self):

dask_geomodeling/raster/elemwise.py

@@ -5,7 +5,7 @@
 
 import numpy as np
 
-from dask_geomodeling.utils import get_dtype_max, get_index, GeoTransform
+from dask_geomodeling.utils import get_dtype_max, get_index, GeoTransform, Extent
 
 from .base import RasterBlock, BaseSingle
 
@@ -156,16 +156,12 @@ def geometry(self):
             return
         if len(geometries) == 1:
             return geometries[0]
-        result = geometries[0]
-        sr = result.GetSpatialReference()
+        extent = Extent.from_geometry(geometries[0])
         for geometry in geometries[1:]:
-            if not geometry.GetSpatialReference().IsSame(sr):
-                geometry = geometry.Clone()
-                geometry.TransformTo(sr)
-            result = result.Intersection(geometry)
-        if result.GetArea() == 0.0:
-            return
-        return result
+            extent = extent.intersection(Extent.from_geometry(geometry))
+            if extent is None:
+                return None
+        return extent.as_geometry()
 
     @property
     def projection(self):

dask_geomodeling/raster/misc.py

@@ -126,14 +126,10 @@ def geometry(self):
         result, mask = [x.geometry for x in self.args]
         if result is None or mask is None:
             return
-        sr = result.GetSpatialReference()
-        if not mask.GetSpatialReference().IsSame(sr):
-            mask = mask.Clone()
-            mask.TransformTo(sr)
-        result = result.Intersection(mask)
-        if result.GetArea() == 0.0:
-            return
-        return result
+        extent = utils.Extent.from_geometry(result).intersection(utils.Extent.from_geometry(mask))
+        if extent is None:
+            return None
+        return extent.as_geometry()
 
     @property
     def period(self):

dask_geomodeling/raster/reduction.py

@@ -2,7 +2,7 @@
 Module containing reduction raster blocks.
 """
 import numpy as np
-from dask_geomodeling.utils import filter_none, get_index
+from dask_geomodeling.utils import filter_none, get_index, Extent
 from dask_geomodeling.utils import parse_percentile_statistic
 from .base import RasterBlock
 from .elemwise import BaseElementwise
@@ -182,14 +182,10 @@ def geometry(self):
             return
         elif len(geometries) == 1:
             return geometries[0]
-        result = geometries[0]
-        sr = result.GetSpatialReference()
+        extent = Extent.from_geometry(geometries[0])
         for geometry in geometries[1:]:
-            if not geometry.GetSpatialReference().IsSame(sr):
-                geometry = geometry.Clone()
-                geometry.TransformTo(sr)
-            result = result.Union(geometry)
-        return result
+            extent = extent.union(Extent.from_geometry(geometry))
+        return extent.as_geometry()
 
 
 def wrap_reduction_function(statistic):

dask_geomodeling/raster/sources.py

@@ -146,14 +146,14 @@ def _get_extent(self):
         if not self.data.size:
             return
         bbox = self.geo_transform.get_bbox((0, 0), self.data.shape[1:])
-        return utils.Extent(bbox, utils.get_sr(self.projection))
+        return utils.Extent(bbox, self.projection)
 
     @property
     def extent(self):
         extent = self._get_extent()
         if extent is None:
             return
-        return extent.transformed(utils.EPSG4326).bbox
+        return extent.transformed("EPSG:4326").bbox
 
     @property
     def geometry(self):
@@ -370,11 +370,11 @@ def _get_extent(self):
         bbox = self.geo_transform.get_bbox(
             (0, 0), (self.gdal_dataset.RasterYSize, self.gdal_dataset.RasterXSize)
         )
-        return utils.Extent(bbox, utils.get_sr(self.projection))
+        return utils.Extent(bbox, self.projection)
 
     @property
     def extent(self):
-        extent_epsg4326 = self._get_extent().transformed(utils.EPSG4326)
+        extent_epsg4326 = self._get_extent().transformed("EPSG:4326")
         return extent_epsg4326.bbox
 
     @property

dask_geomodeling/raster/spatial.py

@@ -8,8 +8,6 @@
 from osgeo import ogr
 
 from dask_geomodeling.utils import (
-    EPSG3857,
-    EPSG4326,
     POLYGON,
     get_sr,
     Extent,
@@ -66,8 +64,8 @@ def expand_request_meters(request, radius_m=1):
     # suffixed by _m, in pixels by _px
     if sr.IsGeographic():
         # expand geographic bbox in EPSG3857
-        extent_geom = Extent(bbox, sr)
-        bbox = extent_geom.transformed(EPSG3857).bbox
+        extent_geom = Extent(bbox, request["projection"])
+        bbox = extent_geom.transformed("EPSG:3857").bbox
     else:
         # most Projected projections are in meters, but to be sure:
         radius_m /= sr.GetLinearUnits()
@@ -102,8 +100,8 @@ def expand_request_meters(request, radius_m=1):
     )
     if sr.IsGeographic():
         # transform back to original projection
-        extent_proj = Extent(new_request["bbox"], EPSG3857)
-        new_request["bbox"] = extent_proj.transformed(sr).bbox
+        extent_proj = Extent(new_request["bbox"], "EPSG:3857")
+        new_request["bbox"] = extent_proj.transformed(request["projection"]).bbox
     new_request["height"] += 2 * margins_px[0]
     new_request["width"] += 2 * margins_px[1]
 
@@ -533,28 +531,21 @@ def extent(self):
         geometry = self.geometry
         if geometry is None:
             return
-        if not geometry.GetSpatialReference().IsSame(EPSG4326):
-            geometry = geometry.Clone()
-            geometry.TransformTo(EPSG4326)
-        x1, x2, y1, y2 = geometry.GetEnvelope()
-        return x1, y1, x2, y2
+        return Extent.from_geometry(geometry).transformed("EPSG:4326").bbox
 
     @property
     def geometry(self):
         """Combined geometry in this block's native projection. """
         store_geometry = self.store.geometry
         if store_geometry is None:
             return
-        sr = get_sr(self.place_projection)
-        if not store_geometry.GetSpatialReference().IsSame(sr):
-            store_geometry = store_geometry.Clone()
-            store_geometry.TransformTo(sr)
-        _x1, _x2, _y1, _y2 = store_geometry.GetEnvelope()
+        extent = Extent.from_geometry(store_geometry).transformed(self.place_projection)
+        _x1, _y1, _x2, _y2 = extent.bbox
         p, q = self.anchor
         P, Q = zip(*self.coordinates)
         x1, x2 = _x1 + min(P) - p, _x2 + max(P) - p
         y1, y2 = _y1 + min(Q) - q, _y2 + max(Q) - q
-        return ogr.CreateGeometryFromWkt(POLYGON.format(x1, y1, x2, y2), sr)
+        return ogr.CreateGeometryFromWkt(POLYGON.format(x1, y1, x2, y2), extent.sr)
 
     def get_sources_and_requests(self, **request):
         if request["mode"] != "vals":
@@ -576,11 +567,7 @@ def get_sources_and_requests(self, **request):
         if extent_geometry is None:
             # no geometry means: no data
             return (({"mode": "null"}, None),)
-        sr = get_sr(request["projection"])
-        if not extent_geometry.GetSpatialReference().IsSame(sr):
-            extent_geometry = extent_geometry.Clone()
-            extent_geometry.TransformTo(sr)
-        xmin, xmax, ymin, ymax = extent_geometry.GetEnvelope()
+        xmin, xmax, ymin, ymax = Extent.from_geometry(extent_geometry).transformed(request["projection"]).bbox
 
         # compute the requested cellsize
         x1, y1, x2, y2 = request["bbox"]

dask_geomodeling/tests/factories.py

@@ -16,7 +16,6 @@
     get_dtype_max,
     Extent,
     get_crs,
-    get_sr,
     get_epsg_or_wkt,
     shapely_transform,
     Dataset,
@@ -114,8 +113,8 @@ def process(args, request):
         bbox = request.get("bbox", (0, 0, width, height))
         projection = request.get("projection", "EPSG:3857")
         if projection != src_projection:
-            extent = Extent(bbox, get_sr(projection))
-            bbox = extent.transformed(get_sr(src_projection)).bbox
+            extent = Extent(bbox, projection)
+            bbox = extent.transformed(src_projection).bbox
         x1, y1, x2, y2 = [int(round(x)) for x in bbox]
 
         if x1 == x2 or y1 == y2:  # point request
@@ -180,7 +179,7 @@ def geo_transform(self):
     def geometry(self):
         if self.extent is None:
             return
-        return Extent(self.extent, get_sr(self.projection)).as_geometry()
+        return Extent(self.extent, self.projection).as_geometry()
 
 
 class MockGeometry(GeometryBlock):

dask_geomodeling/tests/test_geometry.py

@@ -152,8 +152,8 @@ def test_get_data_centroid_mode(self):
         self.assertEqual(0, len(result["features"]))
 
     def test_reproject(self):
-        extent = Extent(self.bbox, get_sr(self.projection))
-        bbox3857 = extent.transformed(get_sr("EPSG:3857")).bbox
+        extent = Extent(self.bbox, self.projection)
+        bbox3857 = extent.transformed("EPSG:3857").bbox
         result = self.source.get_data(geometry=box(*bbox3857), projection="EPSG:3857")
         self.assertEqual("EPSG:3857", result["projection"])
         actual = result["features"].crs
@@ -821,10 +821,10 @@ def test_raster_time_resolution(self):
         req = self.request
         req["time_resolution"] = 3600000
         req["geometry"] = box(0, 0, 10, 10)
-        
+
         # temp_group = GroupTemporal([
         #     MockRaster(timedelta=Timedelta(hours=1)),
-        #     MockRaster(timedelta=Timedelta(minutes=1)), 
+        #     MockRaster(timedelta=Timedelta(minutes=1)),
         #     MockRaster(timedelta=Timedelta(seconds=1))
         # ])
 

dask_geomodeling/tests/test_raster.py

@@ -7,7 +7,6 @@
 from scipy import ndimage
 
 from dask_geomodeling import raster
-from dask_geomodeling.utils import EPSG4326, EPSG3857
 from dask_geomodeling.utils import Extent, get_dtype_max, get_epsg_or_wkt
 from dask_geomodeling.raster import RasterBlock
 from dask_geomodeling.tests.factories import MockRaster, MockGeometry
@@ -1919,8 +1918,8 @@ def test_smooth(self):
         ):  # partial
             request["height"] = bbox[3] - bbox[1]
             request["width"] = bbox[2] - bbox[0]
-            extent = Extent(bbox, EPSG3857)
-            request["bbox"] = extent.transformed(EPSG4326).bbox
+            extent = Extent(bbox, "EPSG:3857")
+            request["bbox"] = extent.transformed("EPSG:4326").bbox
             data = view.get_data(**request)
             _expected = expected[bbox[1] : bbox[3], bbox[0] : bbox[2]]
             assert_allclose(data["values"][0], _expected, atol=peak * 0.0001)

dask_geomodeling/tests/test_raster_misc.py

@@ -44,7 +44,7 @@ def test_clip_attrs_intersects(source, empty_source):
     )
     expected_geometry = source.geometry.Intersection(clipping_mask.geometry)
     assert clip.extent == expected_extent
-    assert clip.geometry.ExportToWkt() == expected_geometry.ExportToWkt()
+    assert clip.geometry.GetEnvelope() == expected_geometry.GetEnvelope()
 
 
 def test_clip_attrs_with_reprojection(source, empty_source):

dask_geomodeling/tests/test_raster_sources.py

@@ -37,15 +37,15 @@ def test_fillvalue(self):
 
     def test_extent(self):
         expected = (
-            utils.Extent((136700, 455795, 136705, 455800), utils.get_sr("EPSG:28992"))
-            .transformed(utils.get_sr("EPSG:4326"))
+            utils.Extent((136700, 455795, 136705, 455800), "EPSG:28992")
+            .transformed("EPSG:4326")
             .bbox
         )
         assert_allclose(self.source.extent, expected, atol=1e-10)
 
     def test_geometry(self):
         expected = (
-            utils.Extent((136700, 455795, 136705, 455800), utils.get_sr("EPSG:28992"))
+            utils.Extent((136700, 455795, 136705, 455800), "EPSG:28992")
             .as_geometry()
             .ExportToWkt()
         )