ENH: offset_curve (#229)

CHANGELOG.rst

@@ -9,6 +9,8 @@ Version 0.9 (unreleased)
 * Release the GIL for ``is_geometry()``, ``is_missing()``, and
   ``is_valid_input()`` (#207)
 * Addition of a ``is_ccw()`` function for GEOS >= 3.7 (#201)
+* Addition of a ``minimum_clearance`` function for GEOS >= 3.6.0 (#223)
+* Addition of a ``offset_curve`` function (#229)
 * Added support for pickling to ``Geometry`` objects (#190)
 * Limited the length of geometry repr to 80 characters (#189)
 * Argument in ``line_interpolate_point`` and ``line_locate_point``

pygeos/constructive.py

@@ -10,6 +10,7 @@
     "BufferJoinStyles",
     "boundary",
     "buffer",
+    "offset_curve",
     "centroid",
     "convex_hull",
     "delaunay_triangles",
@@ -64,6 +65,7 @@ def boundary(geometry, **kwargs):
     """
     return lib.boundary(geometry, **kwargs)
 
+
 @multithreading_enabled
 def buffer(
     geometry,
@@ -166,6 +168,69 @@ def buffer(
         **kwargs
     )
 
+
+@multithreading_enabled
+def offset_curve(
+    geometry,
+    distance,
+    quadsegs=8,
+    join_style="round",
+    mitre_limit=5.0,
+    **kwargs
+):
+    """
+    Returns a (Multi)LineString at a distance from the object
+    on its right or its left side.
+
+    For positive distance the offset will be at the left side of
+    the input line and retain the same direction. For a negative
+    distance it will be at the right side and in the opposite
+    direction.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    distance : float or array_like
+        Specifies the offset distance from the input geometry. Negative
+        for right side offset, positive for left side offset.
+    quadsegs : int
+        Specifies the number of linear segments in a quarter circle in the
+        approximation of circular arcs.
+    join_style : {'round', 'bevel', 'sharp'}
+        Specifies the shape of outside corners. 'round' results in
+        rounded shapes. 'bevel' results in a beveled edge that touches the
+        original vertex. 'mitre' results in a single vertex that is beveled
+        depending on the ``mitre_limit`` parameter.
+    mitre_limit : float
+        Crops of 'mitre'-style joins if the point is displaced from the
+        buffered vertex by more than this limit.
+
+    Examples
+    --------
+    >>> line = Geometry("LINESTRING (0 0, 0 2)")
+    >>> offset_curve(line, 2)
+    <pygeos.Geometry LINESTRING (-2 0, -2 2)>
+    >>> offset_curve(line, -2)
+    <pygeos.Geometry LINESTRING (2 2, 2 0)>
+    """
+    if isinstance(join_style, str):
+        join_style = BufferJoinStyles[join_style.upper()].value
+    if not np.isscalar(quadsegs):
+        raise TypeError("quadsegs only accepts scalar values")
+    if not np.isscalar(join_style):
+        raise TypeError("join_style only accepts scalar values")
+    if not np.isscalar(mitre_limit):
+        raise TypeError("mitre_limit only accepts scalar values")
+    return lib.offset_curve(
+        geometry,
+        distance,
+        np.intc(quadsegs),
+        np.intc(join_style),
+        np.double(mitre_limit),
+        **kwargs
+    )
+
+
 @multithreading_enabled
 def centroid(geometry, **kwargs):
     """Computes the geometric center (center-of-mass) of a geometry.

pygeos/test/test_constructive.py

@@ -4,7 +4,7 @@
 
 from pygeos import Geometry, GEOSException
 
-from .common import point, line_string, all_types, empty
+from .common import point, line_string, all_types, empty, empty_line_string
 
 CONSTRUCTIVE_NO_ARGS = (
     pygeos.boundary,
@@ -18,6 +18,7 @@
 
 CONSTRUCTIVE_FLOAT_ARG = (
     pygeos.buffer,
+    pygeos.offset_curve,
     pygeos.delaunay_triangles,
     pygeos.simplify,
     pygeos.voronoi_polygons,
@@ -35,6 +36,10 @@ def test_no_args_array(geometry, func):
 @pytest.mark.parametrize("geometry", all_types)
 @pytest.mark.parametrize("func", CONSTRUCTIVE_FLOAT_ARG)
 def test_float_arg_array(geometry, func):
+    if func is pygeos.offset_curve and pygeos.get_type_id(geometry) not in [1, 2]:
+        with pytest.raises(GEOSException, match="only accept linestrings"):
+            func([geometry, geometry], 0.0)
+        return
     actual = func([geometry, geometry], 0.0)
     assert actual.shape == (2,)
     assert isinstance(actual[0], Geometry)
@@ -178,3 +183,43 @@ def test_make_valid_1d(geom, expected):
 def test_normalize(geom, expected):
     actual = pygeos.normalize(geom)
     assert actual == expected
+
+
+def test_offset_curve_empty():
+    actual = pygeos.offset_curve(empty_line_string, 2.0)
+    assert pygeos.is_empty(actual)
+
+
+def test_offset_curve_distance_array():
+    # check that kwargs are passed through
+    result = pygeos.offset_curve([line_string, line_string], [-2.0, -3.0])
+    assert result[0] == pygeos.offset_curve(line_string, -2.0)
+    assert result[1] == pygeos.offset_curve(line_string, -3.0)
+
+
+def test_offset_curve_kwargs():
+    # check that kwargs are passed through
+    result1 = pygeos.offset_curve(
+        line_string, -2.0, quadsegs=2, join_style="mitre", mitre_limit=2.0
+    )
+    result2 = pygeos.offset_curve(line_string, -2.0)
+    assert result1 != result2
+
+
+def test_offset_curve_non_scalar_kwargs():
+    msg = "only accepts scalar values"
+    with pytest.raises(TypeError, match=msg):
+        pygeos.offset_curve([line_string, line_string], 1, quadsegs=np.array([8, 9]))
+
+    with pytest.raises(TypeError, match=msg):
+        pygeos.offset_curve(
+            [line_string, line_string], 1, join_style=["round", "bevel"]
+        )
+
+    with pytest.raises(TypeError, match=msg):
+        pygeos.offset_curve([line_string, line_string], 1, mitre_limit=[5.0, 6.0])
+
+
+def test_offset_curve_join_style():
+    with pytest.raises(KeyError):
+        pygeos.offset_curve(line_string, 1.0, join_style="nonsense")

src/ufuncs.c

@@ -166,23 +166,23 @@ static PyUFuncGenericFunction Y_b_funcs[1] = {&Y_b_func};
 /* Define the object -> bool functions (O_b) which do not raise on non-geom objects*/
 static char IsMissing(void* context, PyObject* obj) {
   GEOSGeometry* g = NULL;
-  if (!get_geom((GeometryObject *) obj, &g)) {
+  if (!get_geom((GeometryObject*)obj, &g)) {
     return 0;
   };
   return g == NULL;  // get_geom sets g to NULL for None input
 }
 static void* is_missing_data[1] = {IsMissing};
 static char IsGeometry(void* context, PyObject* obj) {
   GEOSGeometry* g = NULL;
-  if (!get_geom((GeometryObject *) obj, &g)) {
+  if (!get_geom((GeometryObject*)obj, &g)) {
     return 0;
   }
   return g != NULL;
 }
 static void* is_geometry_data[1] = {IsGeometry};
 static char IsValidInput(void* context, PyObject* obj) {
   GEOSGeometry* g = NULL;
-  return get_geom((GeometryObject *) obj, &g);
+  return get_geom((GeometryObject*)obj, &g);
 }
 static void* is_valid_input_data[1] = {IsValidInput};
 typedef char FuncGEOS_O_b(void* context, PyObject* obj);
@@ -1122,6 +1122,79 @@ static void buffer_func(char** args, npy_intp* dimensions, npy_intp* steps, void
 }
 static PyUFuncGenericFunction buffer_funcs[1] = {&buffer_func};
 
+static char offset_curve_dtypes[6] = {NPY_OBJECT, NPY_DOUBLE, NPY_INT,
+                                      NPY_INT,    NPY_DOUBLE, NPY_OBJECT};
+static void offset_curve_func(char** args, npy_intp* dimensions, npy_intp* steps,
+                              void* data) {
+  char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4];
+  npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4];
+  npy_intp n = dimensions[0];
+  npy_intp i;
+  GEOSGeometry** geom_arr;
+  GEOSGeometry* in1 = NULL;
+
+  // Fail if inputs output multiple times on the same place in memory. That would
+  // lead to segfaults as the same GEOSGeometry would be 'owned' by multiple PyObjects.
+  if ((steps[5] == 0) && (dimensions[0] > 1)) {
+    PyErr_Format(PyExc_NotImplementedError,
+                 "Unknown ufunc mode with args[0]=%p, args[5]=%p, steps[0]=%ld, "
+                 "steps[5]=%ld, dimensions[0]=%ld.",
+                 args[0], args[5], steps[0], steps[5], dimensions[0]);
+    return;
+  }
+
+  if ((is3 != 0) | (is4 != 0) | (is5 != 0)) {
+    PyErr_Format(PyExc_ValueError,
+                 "Offset curve function called with non-scalar parameters");
+    return;
+  }
+
+  double width;
+  int quadsegs = *(int*)ip3;
+  int joinStyle = *(int*)ip4;
+  double mitreLimit = *(double*)ip5;
+
+  // allocate a temporary array to store output GEOSGeometry objects
+  geom_arr = malloc(sizeof(void*) * n);
+  if (geom_arr == NULL) {
+    PyErr_SetString(PyExc_MemoryError, "Could not allocate memory");
+    return;
+  }
+
+  GEOS_INIT_THREADS;
+
+  for (i = 0; i < n; i++, ip1 += is1, ip2 += is2) {
+    /* get the geometry: return on error */
+    if (!get_geom(*(GeometryObject**)ip1, &in1)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+
+    width = *(double*)ip2;
+    if ((in1 == NULL) | npy_isnan(width)) {
+      // in case of a missing value: return NULL (None)
+      geom_arr[i] = NULL;
+    } else {
+      geom_arr[i] = GEOSOffsetCurve_r(ctx, in1, width, quadsegs, joinStyle, mitreLimit);
+      if (geom_arr[i] == NULL) {
+        errstate = PGERR_GEOS_EXCEPTION;
+        destroy_geom_arr(ctx, geom_arr, i - 1);
+        break;
+      }
+    }
+  }
+
+  GEOS_FINISH_THREADS;
+
+  // fill the numpy array with PyObjects while holding the GIL
+  if (errstate == PGERR_SUCCESS) {
+    geom_arr_to_npy(geom_arr, args[5], steps[5], dimensions[0]);
+  }
+  free(geom_arr);
+}
+static PyUFuncGenericFunction offset_curve_funcs[1] = {&offset_curve_func};
+
 static char snap_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE, NPY_OBJECT};
 static void snap_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
   GEOSGeometry *in1 = NULL, *in2 = NULL;
@@ -2244,6 +2317,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_YYd_d(hausdorff_distance_densify);
 
   DEFINE_CUSTOM(buffer, 7);
+  DEFINE_CUSTOM(offset_curve, 5);
   DEFINE_CUSTOM(snap, 3);
   DEFINE_CUSTOM(equals_exact, 3);