Solve performance regression with GDAL 3 (#99)

.github/workflows/test.yml

@@ -15,14 +15,10 @@ jobs:
       fail-fast: false
       matrix:
         include:
-          - os: ubuntu-18.04
-            python: 3.6
-            numpy: "==1.14.*"
-            pins: "pygdal==2.2.3.* scipy==1.1.* dask[delayed]==0.20.* pandas==0.23.* geopandas==0.7.*"
-          - os: ubuntu-18.04
+          - os: ubuntu-20.04
             python: 3.7
             numpy: "==1.16.*"
-            pins: "pygdal==2.2.3.* scipy==1.3.* dask[delayed]==1.* pandas==0.25.* geopandas==0.7.*"
+            pins: "pygdal==3.0.4.* scipy==1.3.* dask[delayed]==1.* pandas==0.25.* geopandas==0.7.*"
           - os: ubuntu-20.04
             python: 3.8
             numpy: "==1.18.*"
@@ -40,7 +36,7 @@ jobs:
             numpy: "==1.23.*"
             pins: "pygdal==3.4.1.* scipy==1.9.* dask[delayed]==2021.7.* pandas==1.4.* geopandas==0.11.*"
           - os: ubuntu-22.04
-            python: "3.10"
+            python: "3.11"
             display_name: "latest"
             pins: "pygdal==3.4.1.*"
 

CHANGES.rst

@@ -4,7 +4,10 @@ Changelog of dask-geomodeling
 2.3.11 (unreleased)
 -------------------
 
-- Nothing changed yet.
+- Perform geometry transformations via pyproj instead of GDAL SWIG bindings to circumvent
+  a performance degradation in GDAL >=3 (PROJ >=6).
+
+- Drop support for Python 3.6 and GDAL 2.
 
 
 2.3.10 (2022-08-22)

dask_geomodeling/tests/test_utils.py

@@ -16,6 +16,8 @@
 
 from dask_geomodeling import utils
 
+from osgeo import ogr
+
 
 class TestUtils(unittest.TestCase):
     def test_get_index(self):
@@ -190,14 +192,21 @@ def test_get_crs(self):
     def test_shapely_transform(self):
         src_srs = "EPSG:28992"
         dst_srs = "EPSG:4326"
-        box28992 = geometry.box(100000, 400000, 200000, 500000)
+        box28992 = geometry.box(100000, 400000, 101000, 401000)
         box4326 = utils.shapely_transform(box28992, src_srs=src_srs, dst_srs=dst_srs)
-        assert_almost_equal((5.4, 52.0), list(box4326.centroid.coords)[0], decimal=1)
+        assert_almost_equal((4.608024, 51.586315), box4326.exterior.coords[0], decimal=6)
 
     def test_shapely_transform_invalid(self):
         src_srs = "EPSG:4326"
         dst_srs = "EPSG:28992"
-        box4326 = geometry.box(100000, 400000, 200000, 500000)
+        box4326 = geometry.box(100000, 400000, 101000, 401000)
+        with pytest.raises(utils.TransformException):
+            utils.shapely_transform(box4326, src_srs=src_srs, dst_srs=dst_srs)
+
+    def test_shapely_transform_srs(self):
+        src_srs = "EPSG:0"
+        dst_srs = "EPSG:28992"
+        box4326 = geometry.box(100000, 400000, 101000, 401000)
         with pytest.raises(utils.TransformException):
             utils.shapely_transform(box4326, src_srs=src_srs, dst_srs=dst_srs)
 

dask_geomodeling/utils.py

@@ -3,7 +3,7 @@
 import pytz
 import os
 import warnings
-from functools import lru_cache
+from functools import lru_cache, partial
 from itertools import repeat
 
 from math import floor, log10
@@ -15,6 +15,7 @@
 
 from osgeo import gdal, ogr, osr, gdal_array
 from shapely.geometry import box, Point
+from shapely.ops import transform as _shapely_transform
 
 try:  # shapely 2.*
     from shapely import from_wkt, from_wkb, GEOSException
@@ -24,7 +25,8 @@
     from shapely.wkb import loads as from_wkb
 
 
-from pyproj import CRS
+from pyproj import CRS, Transformer
+from pyproj.exceptions import ProjError
 
 import fiona
 import fiona.crs
@@ -403,18 +405,12 @@ class TransformException(Exception):
     pass
 
 
-def wkb_transform(wkb, src_sr, dst_sr):
-    """
-    Return a shapely geometry transformed from src_sr to dst_sr.
-
-    :param wkb: wkb bytes
-    :param src_sr: source osr SpatialReference
-    :param dst_sr: destination osr SpatialReference
-    """
-    result = ogr.CreateGeometryFromWkb(wkb, src_sr)
-    result.TransformTo(dst_sr)
-    wkb = result.ExportToWkb()
-    return bytes(wkb) if not isinstance(wkb, bytes) else wkb
+@lru_cache(maxsize=100)
+def get_transform_func(src_srs: str, dst_srs: str):
+    transformer = Transformer.from_crs(
+        CRS(src_srs), CRS(dst_srs), always_xy=True
+    )
+    return partial(transformer.transform, errcheck=True)
 
 
 def shapely_transform(geometry, src_srs, dst_srs):
@@ -425,27 +421,18 @@ def shapely_transform(geometry, src_srs, dst_srs):
     :param src_srs: source projection string
     :param dst_srs: destination projection string
 
-    Note that we do not use geopandas for the transformation, because this is
-    much slower than OGR.
+    Note that we separately construct (and cache) the Transformer instance,
+    because of large overhead in constructing the instance since PROJ v6. 
     """
-    transform_kwargs = {
-        "wkb": geometry.wkb,
-        "src_sr": get_sr(src_srs),
-        "dst_sr": get_sr(dst_srs),
-    }
     try:
-        output_wkb = wkb_transform(**transform_kwargs)
-    except RuntimeError:
-        # include the geometry in the error message, truncated to 64 chars
-        wkt = geometry.wkt
-        if len(wkt) > 64:
-            wkt = wkt[:61] + "..."
+        func = get_transform_func(src_srs, dst_srs)
+        return _shapely_transform(func, geometry)
+    except ProjError:
         raise TransformException(
             "An error occured while transforming {} from {} to {}.".format(
-                wkt, src_srs, dst_srs
+                geometry.wkt, src_srs, dst_srs
             )
         )
-    return from_wkb(output_wkb)
 
 
 def shapely_from_wkt(wkt):

setup.py

@@ -52,7 +52,7 @@
     zip_safe=False,
     install_requires=install_requires,
     tests_require=tests_require,
-    python_requires='>=3.6',
+    python_requires='>=3.7',
     extras_require={"test": tests_require, "cityhash": ["cityhash"]},
     entry_points={"console_scripts": []},
 )