Python Bindings for Pairwise Linestring Distance

cpp/include/cuspatial/distances/linestring_distance.hpp

@@ -83,8 +83,8 @@ namespace cuspatial {
  * @param linestring1_points_x x-components of points in the first linestring of each pair.
  * @param linestring1_points_y y-component of points in the first linestring of each pair.
  * @param linestring2_offsets Indices of the first point of the second linestring of each pair.
- * @param linestring2_points_x x-component of points in the first linestring of each pair.
- * @param linestring2_points_y y-component of points in the first linestring of each pair.
+ * @param linestring2_points_x x-component of points in the second linestring of each pair.
+ * @param linestring2_points_y y-component of points in the second linestring of each pair.
  * @param mr Device memory resource used to allocate the returned column's device memory.
  * @return A column of shortest distances between each pair of linestrings.
  *

python/cuspatial/cuspatial/__init__.py

@@ -7,6 +7,7 @@
     point_in_polygon,
     polygon_bounding_boxes,
     polyline_bounding_boxes,
+    pairwise_linestring_distance,
 )
 from .core.indexing import quadtree_on_points
 from .core.interpolate import CubicSpline

python/cuspatial/cuspatial/_lib/cpp/linestring_distance.pxd

@@ -0,0 +1,20 @@
+# Copyright (c) 2022, NVIDIA CORPORATION.
+
+from libcpp.memory cimport unique_ptr
+
+from cudf._lib.column cimport Column
+from cudf._lib.cpp.column.column cimport column
+from cudf._lib.cpp.column.column_view cimport column_view
+from cudf._lib.cpp.types cimport size_type
+
+
+cdef extern from "cuspatial/distances/linestring_distance.hpp" \
+        namespace "cuspatial" nogil:
+    cdef unique_ptr[column] pairwise_linestring_distance(
+        const column_view linestring1_offsets,
+        const column_view linestring1_points_x,
+        const column_view linestring1_points_y,
+        const column_view linestring2_offsets,
+        const column_view linestring2_points_x,
+        const column_view linestring2_points_y
+    ) except +

python/cuspatial/cuspatial/_lib/linestring_distance.pyx

@@ -0,0 +1,39 @@
+from libcpp.memory cimport unique_ptr
+from libcpp.utility cimport move
+
+from cudf._lib.column cimport Column
+from cudf._lib.cpp.column.column cimport column
+from cudf._lib.cpp.column.column_view cimport column_view
+
+from cuspatial._lib.cpp.linestring_distance cimport (
+    pairwise_linestring_distance as cpp_pairwise_linestring_distance,
+)
+
+
+def pairwise_linestring_distance(
+    Column linestring1_offsets,
+    Column linestring1_points_x,
+    Column linestring1_points_y,
+    Column linestring2_offsets,
+    Column linestring2_points_x,
+    Column linestring2_points_y
+):
+    cdef column_view linestring1_offsets_view = linestring1_offsets.view()
+    cdef column_view linestring1_points_x_view = linestring1_points_x.view()
+    cdef column_view linestring1_points_y_view = linestring1_points_y.view()
+    cdef column_view linestring2_offsets_view = linestring2_offsets.view()
+    cdef column_view linestring2_points_x_view = linestring2_points_x.view()
+    cdef column_view linestring2_points_y_view = linestring2_points_y.view()
+
+    cdef unique_ptr[column] c_result
+    with nogil:
+        c_result = move(cpp_pairwise_linestring_distance(
+            linestring1_offsets_view,
+            linestring1_points_x_view,
+            linestring1_points_y_view,
+            linestring2_offsets_view,
+            linestring2_points_x_view,
+            linestring2_points_y_view
+        ))
+
+    return Column.from_unique_ptr(move(c_result))

python/cuspatial/cuspatial/core/gis.py

@@ -1,11 +1,14 @@
 # Copyright (c) 2019-2020, NVIDIA CORPORATION.
 
-from cudf import DataFrame
+from cudf import DataFrame, Series
 from cudf.core.column import as_column
 
 from cuspatial._lib.hausdorff import (
     directed_hausdorff_distance as cpp_directed_hausdorff_distance,
 )
+from cuspatial._lib.linestring_distance import (
+    pairwise_linestring_distance as cpp_pairwise_linestring_distance,
+)
 from cuspatial._lib.point_in_polygon import (
     point_in_polygon as cpp_point_in_polygon,
 )
@@ -335,3 +338,102 @@ def polyline_bounding_boxes(poly_offsets, xs, ys, expansion_radius):
     return DataFrame._from_data(
         *cpp_polyline_bounding_boxes(poly_offsets, xs, ys, expansion_radius)
     )
+
+
+def pairwise_linestring_distance(offsets1, xs1, ys1, offsets2, xs2, ys2):
+    """Compute shortest distance between pairs of linestrings (a.k.a. polylines)
+
+    The shortest distance between two linestrings is defined as the shortest
+    distance between all pairs of segments of the two linestrings. If any of
+    the segments intersect, the distance is 0.
+
+    Parameters
+    ----------
+    offsets1
+        Indices of the first point of the first linestring of each pair.
+    xs1
+        x-components of points in the first linestring of each pair.
+    ys1
+        y-component of points in the first linestring of each pair.
+    offsets2
+        Indices of the first point of the second linestring of each pair.
+    xs2
+        x-component of points in the second linestring of each pair.
+    ys2
+        y-component of points in the second linestring of each pair.
+
+    Returns
+    -------
+    distance : cudf.Series
+        the distance between each pair of linestrings
+
+    Examples
+    --------
+    The following example contains 4 pairs of linestrings.
+
+    First pair::
+
+        (0, 1) -> (1, 0) -> (-1, 0)
+        (1, 1) -> (2, 1) -> (2, 0) -> (3, 0)
+
+            |
+            *   #####
+            | *     #
+        ----O---*---#####
+            | *
+            *
+            |
+
+    The shortest distance between the two linestrings is the distance
+    from point ``(1, 1)`` to segment ``(0, 1) -> (1, 0)``, which is
+    ``sqrt(2)/2``.
+
+    Second pair::
+
+        (0, 0) -> (0, 1)
+        (1, 0) -> (1, 1) -> (1, 2)
+
+
+    These linestrings are parallel. Their distance is 1 (point
+    ``(0, 0)`` to point ``(1, 0)``).
+
+    Third pair::
+
+        (0, 0) -> (2, 2) -> (-2, 0)
+        (2, 0) -> (0, 2)
+
+
+    These linestrings intersect, so their distance is 0.
+
+    Forth pair::
+
+        (2, 2) -> (-2, -2)
+        (1, 1) -> (5, 5) -> (10, 0)
+
+
+    These linestrings contain colinear and overlapping sections, so
+    their distance is 0.
+
+    The input of above example is::
+
+        linestring1_offsets:  {0, 3, 5, 8}
+        linestring1_points_x: {0, 1, -1, 0, 0, 0, 2, -2, 2, -2}
+        linestring1_points_y: {1, 0, 0, 0, 1, 0, 2, 0, 2, -2}
+        linestring2_offsets:  {0, 4, 7, 9}
+        linestring2_points_x: {1, 2, 2, 3, 1, 1, 1, 2, 0, 1, 5, 10}
+        linestring2_points_y: {1, 1, 0, 0, 0, 1, 2, 0, 2, 1, 5, 0}
+
+        Result: {sqrt(2.0)/2, 1, 0, 0}
+    """
+    xs1, ys1, xs2, ys2 = normalize_point_columns(
+        as_column(xs1), as_column(ys1), as_column(xs2), as_column(ys2)
+    )
+    offsets1 = as_column(offsets1, dtype="int32")
+    offsets2 = as_column(offsets2, dtype="int32")
+    return Series._from_data(
+        {
+            None: cpp_pairwise_linestring_distance(
+                offsets1, xs1, ys1, offsets2, xs2, ys2
+            )
+        }
+    )

python/cuspatial/cuspatial/tests/test_pairwise_linestring_distance.py

@@ -0,0 +1,152 @@
+import cupy as cp
+import pandas as pd
+import shapely
+
+import cudf
+from cudf.testing._utils import assert_eq
+
+import cuspatial
+
+
+def shapely_pairwise_linestring_distance(data1, data2, offset1, offset2):
+    """Compute pairwise linestring distances with shapely."""
+
+    def make_linestring(group):
+        return shapely.geometry.LineString([*zip(group["x"], group["y"])])
+
+    ridx1 = pd.RangeIndex(len(data1))
+    ridx2 = pd.RangeIndex(len(data2))
+    groupid1 = ridx1.map(lambda i: offset1.searchsorted(i, side="right"))
+    groupid2 = ridx2.map(lambda i: offset2.searchsorted(i, side="right"))
+
+    data1["gid"] = groupid1
+    data2["gid"] = groupid2
+
+    linestrings1 = data1.groupby("gid").apply(make_linestring)
+    linestrings2 = data2.groupby("gid").apply(make_linestring)
+
+    linestring_pairs = pd.DataFrame({"s1": linestrings1, "s2": linestrings2})
+    distances = linestring_pairs.apply(
+        lambda row: row["s1"].distance(row["s2"]), axis=1
+    )
+
+    return distances.reset_index(drop=True)
+
+
+def test_zero_pair():
+    data1 = cudf.DataFrame(
+        {
+            "x": [],
+            "y": [],
+        }
+    )
+    data2 = cudf.DataFrame(
+        {
+            "x": [],
+            "y": [],
+        }
+    )
+    offset1 = cudf.Series([], dtype="int32")
+    offset2 = cudf.Series([], dtype="int32")
+
+    got = cuspatial.pairwise_linestring_distance(
+        data1["x"], data2["y"], offset1, data2["x"], data2["y"], offset2
+    )
+    expected = cudf.Series([], dtype="float64")
+
+    assert_eq(got, expected)
+
+
+def test_one_pair():
+    data1 = cudf.DataFrame(
+        {
+            "x": [0.0, 1.0],
+            "y": [0.0, 1.0],
+        }
+    )
+    data2 = cudf.DataFrame(
+        {
+            "x": [2.0, 3.0],
+            "y": [2.0, 3.0],
+        }
+    )
+    offset1 = cudf.Series([0], dtype="int32")
+    offset2 = cudf.Series([0], dtype="int32")
+
+    got = cuspatial.pairwise_linestring_distance(
+        data1["x"], data2["y"], offset1, data2["x"], data2["y"], offset2
+    )
+    expected = shapely_pairwise_linestring_distance(
+        data1.to_pandas(),
+        data2.to_pandas(),
+        offset1.to_pandas(),
+        offset2.to_pandas(),
+    )
+
+    assert_eq(got, expected)
+
+
+def test_two_pairs():
+    data1 = cudf.DataFrame(
+        {
+            "x": [0.0, 1.0, 5.0, 7.0, 8.0],
+            "y": [0.0, 1.0, 10.2, 11.4, 12.8],
+        }
+    )
+    data2 = cudf.DataFrame(
+        {
+            "x": [2.0, 3.0, -8.0, -10.0, -13.0, -3.0],
+            "y": [2.0, 3.0, -8.0, -5.0, -15.0, -6.0],
+        }
+    )
+    offset1 = cudf.Series([0, 3], dtype="int32")
+    offset2 = cudf.Series([0, 2], dtype="int32")
+
+    got = cuspatial.pairwise_linestring_distance(
+        data1["x"], data1["y"], offset1, data2["x"], data2["y"], offset2
+    )
+    expected = shapely_pairwise_linestring_distance(
+        data1.to_pandas(),
+        data2.to_pandas(),
+        offset1.to_pandas(),
+        offset2.to_pandas(),
+    )
+
+    assert_eq(got, expected)
+
+
+def test_100_randomized_input():
+    rng = cp.random.RandomState(0)
+
+    max_linestring_points = 10
+    size = 100
+
+    offset1 = rng.randint(2, max_linestring_points, size=(size,))
+    offset2 = rng.randint(2, max_linestring_points, size=(size,))
+
+    offset1 = cp.cumsum(offset1)
+    offset2 = cp.cumsum(offset2)
+
+    num_points_1 = int(offset1[-1])
+    num_points_2 = int(offset2[-1])
+
+    offset1 = cp.concatenate((cp.zeros((1,)), offset1[:-1]))
+    offset2 = cp.concatenate((cp.zeros((1,)), offset2[:-1]))
+
+    points1_x = rng.uniform(-1, 1, (num_points_1,))
+    points1_y = rng.uniform(-1, 1, (num_points_1,))
+
+    points2_x = rng.uniform(0.5, 2.5, (num_points_2,))
+    points2_y = rng.uniform(0.5, 2.5, (num_points_2,))
+
+    got = cuspatial.pairwise_linestring_distance(
+        points1_x, points1_y, offset1, points2_x, points2_y, offset2
+    )
+    expected = shapely_pairwise_linestring_distance(
+        pd.DataFrame({"x": points1_x.get(), "y": points1_y.get()}),
+        pd.DataFrame({"x": points2_x.get(), "y": points2_y.get()}),
+        pd.Series(offset1.get()),
+        pd.Series(offset2.get()),
+    )
+
+    assert_eq(got, expected)