ENH: Add densify ufunc (#299)

CHANGELOG.rst

@@ -11,14 +11,16 @@ Version 0.10 (unreleased)
 **API Changes**
 
 * STRtree default leaf size is now 10 instead of 5, for somewhat better performance
-  under normal conditions.
+  under normal conditions (#286)
 
 **Added GEOS functions**
 
-* Addition of a ``contains_properly`` function (#267).
-* Addition of a ``polygonize`` function (#275).
+* Addition of a ``contains_properly`` function (#267)
+* Addition of a ``polygonize`` function (#275)
+* Addition of a ``segmentize`` function for GEOS >= 3.10 (#299)
 
 **Bug fixes**
+
 * Fixed portability issue for ARM architecture (#293)
 
 **Acknowledgments**

pygeos/constructive.py

@@ -15,6 +15,7 @@
     "clip_by_rect",
     "convex_hull",
     "delaunay_triangles",
+    "segmentize",
     "envelope",
     "extract_unique_points",
     "build_area",
@@ -566,6 +567,39 @@ def reverse(geometry, **kwargs):
     return lib.reverse(geometry, **kwargs)
 
 
+@requires_geos("3.10.0")
+@multithreading_enabled
+def segmentize(geometry, tolerance, **kwargs):
+    """Adds vertices to line segments based on tolerance.
+
+    Additional vertices will be added to every line segment in an input geometry
+    so that segments are no greater than tolerance.  New vertices will evenly
+    subdivide each segment.
+
+    Only linear components of input geometries are densified; other geometries
+    are returned unmodified.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    tolerance : float or array_like
+        Additional vertices will be added so that all line segments are no
+        greater than this value.  Must be greater than 0.
+
+    Examples
+    --------
+    >>> line = Geometry("LINESTRING (0 0, 0 10)")
+    >>> segmentize(line, tolerance=5)
+    <pygeos.Geometry LINESTRING (0 0, 0 5, 0 10)>
+    >>> poly = Geometry("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))")
+    >>> segmentize(poly, tolerance=5)
+    <pygeos.Geometry POLYGON ((0 0, 5 0, 10 0, 10 5, 10 10, 5 10, 0 10, 0 5, 0 0))>
+    >>> segmentize(None, tolerance=5) is None
+    True
+    """
+    return lib.segmentize(geometry, tolerance, **kwargs)
+
+
 @multithreading_enabled
 def simplify(geometry, tolerance, preserve_topology=False, **kwargs):
     """Returns a simplified version of an input geometry using the

pygeos/test/test_constructive.py

@@ -340,13 +340,13 @@ def test_clip_by_rect(geom, expected):
         (
             "POLYGON ((0 0, 0 30, 30 30, 30 0, 0 0), (10 10, 20 10, 20 20, 10 20, 10 10))",
             (10, 10, 20, 20),
-            "GEOMETRYCOLLECTION EMPTY"
+            "GEOMETRYCOLLECTION EMPTY",
         ),
         # Polygon hole (CW) fully on rectangle boundary"""
         (
             "POLYGON ((0 0, 0 30, 30 30, 30 0, 0 0), (10 10, 10 20, 20 20, 20 10, 10 10))",
             (10, 10, 20, 20),
-            "GEOMETRYCOLLECTION EMPTY"
+            "GEOMETRYCOLLECTION EMPTY",
         ),
         # Polygon fully within rectangle"""
         (
@@ -359,7 +359,7 @@ def test_clip_by_rect(geom, expected):
             "POLYGON ((0 0, 0 30, 30 30, 30 0, 0 0), (10 10, 20 10, 20 20, 10 20, 10 10))",
             (5, 5, 15, 15),
             "POLYGON ((5 5, 5 15, 10 15, 10 10, 15 10, 15 5, 5 5))",
-        )
+        ),
     ],
 )
 def test_clip_by_rect_polygon(geom, rect, expected):
@@ -468,3 +468,114 @@ def test_polygonize_missing():
     # set of geometries that is all missing
     result = pygeos.polygonize([None, None])
     assert result == pygeos.Geometry("GEOMETRYCOLLECTION EMPTY")
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize("geometry", all_types)
+@pytest.mark.parametrize("tolerance", [-1, 0])
+def test_segmentize_invalid_tolerance(geometry, tolerance):
+    with pytest.raises(GEOSException, match="IllegalArgumentException"):
+        pygeos.segmentize(geometry, tolerance=tolerance)
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize("geometry", all_types)
+def test_segmentize_tolerance_nan(geometry):
+    actual = pygeos.segmentize(geometry, tolerance=np.nan)
+    assert actual == None
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize(
+    "geometry",
+    [
+        empty,
+        empty_point,
+        empty_line_string,
+        empty_polygon,
+    ],
+)
+def test_segmentize_empty(geometry):
+    actual = pygeos.segmentize(geometry, tolerance=5)
+    assert pygeos.equals(actual, geometry).all()
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize("geometry", [point, point_z, multi_point])
+def test_segmentize_no_change(geometry):
+    actual = pygeos.segmentize(geometry, tolerance=5)
+    assert pygeos.equals(actual, geometry).all()
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+def test_segmentize_none():
+    assert pygeos.segmentize(None, tolerance=5) is None
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize(
+    "geometry,tolerance, expected",
+    [
+        # tolerance greater than max edge length, no change
+        (
+            pygeos.Geometry("LINESTRING (0 0, 0 10)"),
+            20,
+            pygeos.Geometry("LINESTRING (0 0, 0 10)"),
+        ),
+        (
+            pygeos.Geometry("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))"),
+            20,
+            pygeos.Geometry("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))"),
+        ),
+        # tolerance causes one vertex per segment
+        (
+            pygeos.Geometry("LINESTRING (0 0, 0 10)"),
+            5,
+            pygeos.Geometry("LINESTRING (0 0, 0 5, 0 10)"),
+        ),
+        (
+            Geometry("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))"),
+            5,
+            pygeos.Geometry(
+                "POLYGON ((0 0, 5 0, 10 0, 10 5, 10 10, 5 10, 0 10, 0 5, 0 0))"
+            ),
+        ),
+        # ensure input arrays are broadcast correctly
+        (
+            [
+                pygeos.Geometry("LINESTRING (0 0, 0 10)"),
+                pygeos.Geometry("LINESTRING (0 0, 0 2)"),
+            ],
+            5,
+            [
+                pygeos.Geometry("LINESTRING (0 0, 0 5, 0 10)"),
+                pygeos.Geometry("LINESTRING (0 0, 0 2)"),
+            ],
+        ),
+        (
+            [
+                pygeos.Geometry("LINESTRING (0 0, 0 10)"),
+                pygeos.Geometry("LINESTRING (0 0, 0 2)"),
+            ],
+            [5],
+            [
+                pygeos.Geometry("LINESTRING (0 0, 0 5, 0 10)"),
+                pygeos.Geometry("LINESTRING (0 0, 0 2)"),
+            ],
+        ),
+        (
+            [
+                pygeos.Geometry("LINESTRING (0 0, 0 10)"),
+                pygeos.Geometry("LINESTRING (0 0, 0 2)"),
+            ],
+            [5, 1.5],
+            [
+                pygeos.Geometry("LINESTRING (0 0, 0 5, 0 10)"),
+                pygeos.Geometry("LINESTRING (0 0, 0 1, 0 2)"),
+            ],
+        ),
+    ],
+)
+def test_segmentize(geometry, tolerance, expected):
+    actual = pygeos.segmentize(geometry, tolerance)
+    assert pygeos.equals(actual, geometry).all()

src/geos.h

@@ -124,6 +124,7 @@ enum {
 #define GEOS_SINCE_3_7_0 ((GEOS_VERSION_MAJOR >= 3) && (GEOS_VERSION_MINOR >= 7))
 #define GEOS_SINCE_3_8_0 ((GEOS_VERSION_MAJOR >= 3) && (GEOS_VERSION_MINOR >= 8))
 #define GEOS_SINCE_3_9_0 ((GEOS_VERSION_MAJOR >= 3) && (GEOS_VERSION_MINOR >= 9))
+#define GEOS_SINCE_3_10_0 ((GEOS_VERSION_MAJOR >= 3) && (GEOS_VERSION_MINOR >= 10))
 
 extern PyObject* geos_exception[1];
 

src/ufuncs.c

@@ -547,6 +547,10 @@ static void* simplify_preserve_topology_data[1] = {GEOSTopologyPreserveSimplify_
 static void* unary_union_prec_data[1] = {GEOSUnaryUnionPrec_r};
 #endif
 
+#if GEOS_SINCE_3_10_0
+static void* segmentize_data[1] = {GEOSDensify_r};
+#endif
+
 typedef void* FuncGEOS_Yd_Y(void* context, void* a, double b);
 static char Yd_Y_dtypes[3] = {NPY_OBJECT, NPY_DOUBLE, NPY_OBJECT};
 static void Yd_Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
@@ -1189,7 +1193,6 @@ static void YYd_d_func(char** args, npy_intp* dimensions, npy_intp* steps, void*
 }
 static PyUFuncGenericFunction YYd_d_funcs[1] = {&YYd_d_func};
 
-
 #if GEOS_SINCE_3_9_0
 
 /* Define the geom, geom, double -> geom functions (YYd_Y) */
@@ -1244,7 +1247,6 @@ static void YYd_Y_func(char** args, npy_intp* dimensions, npy_intp* steps, void*
 static PyUFuncGenericFunction YYd_Y_funcs[1] = {&YYd_Y_func};
 #endif
 
-
 /* Define functions with unique call signatures */
 
 static void* null_data[1] = {NULL};
@@ -2789,6 +2791,10 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_Yd_Y(unary_union_prec);
 #endif
 
+#if GEOS_SINCE_3_10_0
+  DEFINE_Yd_Y(segmentize);
+#endif
+
   Py_DECREF(ufunc);
   return 0;
 }