Python API for point-linestring nearest points

docs/source/api_docs/spatial.rst

@@ -24,6 +24,11 @@ Measurement Functions
 .. autofunction:: cuspatial.haversine_distance
 .. autofunction:: cuspatial.pairwise_linestring_distance
 
+Nearest Points Function
++++++++++++++++++++++++
+
+.. autofunction:: cuspatial.pairwise_point_linestring_nearest_points
+
 Bounding Boxes
 ++++++++++++++
 

docs/source/api_docs/trajectory.rst

@@ -7,5 +7,4 @@ Functions for identifying and grouping trajectories from point data.
 
 .. autofunction:: cuspatial.derive_trajectories
 .. autofunction:: cuspatial.trajectory_distances_and_speeds
-.. autofunction:: cuspatial.directed_hausdorff_distance
 .. autofunction:: cuspatial.trajectory_bounding_boxes

python/cuspatial/cuspatial/__init__.py

@@ -7,6 +7,7 @@
     lonlat_to_cartesian,
     pairwise_linestring_distance,
     pairwise_point_linestring_distance,
+    pairwise_point_linestring_nearest_points,
     polygon_bounding_boxes,
     polyline_bounding_boxes,
     point_in_polygon,

python/cuspatial/cuspatial/_lib/CMakeLists.txt

@@ -15,6 +15,7 @@
 set(cython_sources 
     hausdorff.pyx
     interpolate.pyx
+    nearest_points.pyx
     point_in_polygon.pyx
     polygon_bounding_boxes.pyx
     polyline_bounding_boxes.pyx

python/cuspatial/cuspatial/_lib/cpp/point_linestring_nearest_points.pxd

@@ -0,0 +1,26 @@
+# Copyright (c) 2022, NVIDIA CORPORATION.
+
+from libcpp.memory cimport unique_ptr
+
+from cudf._lib.column cimport column, column_view
+
+from cuspatial._lib.cpp.optional cimport optional
+
+
+cdef extern from "cuspatial/point_linestring_nearest_points.hpp" \
+        namespace "cuspatial" nogil:
+
+    cdef struct point_linestring_nearest_points_result:
+        optional[unique_ptr[column]] nearest_point_geometry_id
+        optional[unique_ptr[column]] nearest_linestring_geometry_id
+        unique_ptr[column] nearest_segment_id
+        unique_ptr[column] nearest_point_on_linestring_xy
+
+    cdef point_linestring_nearest_points_result \
+        pairwise_point_linestring_nearest_points(
+            const optional[column_view] multipoint_geometry_offsets,
+            const column_view points_xy,
+            const optional[column_view] multilinestring_geometry_offsets,
+            const column_view linestring_part_offsets,
+            const column_view linestring_points_xy,
+        ) except +

python/cuspatial/cuspatial/_lib/nearest_points.pyx

@@ -0,0 +1,78 @@
+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.optional cimport nullopt, optional
+from cuspatial._lib.cpp.point_linestring_nearest_points cimport (
+    pairwise_point_linestring_nearest_points as c_func,
+    point_linestring_nearest_points_result,
+)
+
+
+def pairwise_point_linestring_nearest_points(
+    Column points_xy,
+    Column linestring_part_offsets,
+    Column linestring_points_xy,
+    multipoint_geometry_offset=None,
+    multilinestring_geometry_offset=None,
+):
+    cdef Column multipoint_geometry_offset_column
+    cdef Column multilinestring_geometry_offset_column
+    cdef optional[column_view] c_multipoint_geometry_offset
+    cdef optional[column_view] c_multilinestring_geometry_offset
+
+    if multipoint_geometry_offset is not None:
+        multipoint_geometry_offset_column = multipoint_geometry_offset
+        c_multipoint_geometry_offset = multipoint_geometry_offset_column.view()
+    else:
+        c_multipoint_geometry_offset = nullopt
+
+    if multilinestring_geometry_offset is not None:
+        multilinestring_geometry_offset_column = (
+            multilinestring_geometry_offset
+        )
+        c_multilinestring_geometry_offset = (
+            multilinestring_geometry_offset_column.view()
+        )
+    else:
+        c_multilinestring_geometry_offset = nullopt
+
+    cdef column_view c_points_xy = points_xy.view()
+    cdef column_view c_linestring_offsets = linestring_part_offsets.view()
+    cdef column_view c_linestring_points_xy = linestring_points_xy.view()
+    cdef point_linestring_nearest_points_result c_result
+
+    with nogil:
+        c_result = move(c_func(
+            c_multipoint_geometry_offset,
+            c_points_xy,
+            c_multilinestring_geometry_offset,
+            c_linestring_offsets,
+            c_linestring_points_xy,
+        ))
+
+    if multipoint_geometry_offset is not None:
+        point_geometry_id = Column.from_unique_ptr(
+            move(c_result.nearest_point_geometry_id.value()))
+    else:
+        point_geometry_id = None
+
+    if multilinestring_geometry_offset is not None:
+        linestring_geometry_id = Column.from_unique_ptr(
+            move(c_result.nearest_linestring_geometry_id.value()))
+    else:
+        linestring_geometry_id = None
+
+    segment_id = Column.from_unique_ptr(move(c_result.nearest_segment_id))
+    point_on_linestring_xy = Column.from_unique_ptr(
+        move(c_result.nearest_point_on_linestring_xy))
+
+    return (
+        point_geometry_id,
+        linestring_geometry_id,
+        segment_id,
+        point_on_linestring_xy
+    )

python/cuspatial/cuspatial/core/_column/geocolumn.py

@@ -5,9 +5,9 @@
 import pyarrow as pa
 
 import cudf
-from cudf.core.column import ColumnBase, as_column
+from cudf.core.column import ColumnBase, as_column, build_list_column
 
-from cuspatial.core._column.geometa import GeoMeta
+from cuspatial.core._column.geometa import Feature_Enum, GeoMeta
 
 T = TypeVar("T", bound="GeoColumn")
 
@@ -181,3 +181,40 @@ def copy(self, deep=True):
             self.data.copy(),
         )
         return result
+
+    @classmethod
+    def _from_points_xy(cls, points_xy: ColumnBase):
+        """
+        Create a GeoColumn of only single points from a cudf Series with
+        interleaved xy coordinates.
+        """
+        if len(points_xy) % 2 != 0:
+            raise ValueError("points_xy must have an even number of elements")
+
+        num_points = len(points_xy) // 2
+        meta = GeoMeta(
+            {
+                "input_types": as_column(
+                    cp.full(
+                        num_points, Feature_Enum.POINT.value, dtype=cp.int8
+                    )
+                ),
+                "union_offsets": as_column(
+                    cp.arange(num_points, dtype=cp.int32)
+                ),
+            }
+        )
+
+        indices = as_column(cp.arange(0, num_points * 2 + 1, 2), dtype="int32")
+        point_col = build_list_column(
+            indices=indices, elements=points_xy, size=num_points
+        )
+        return cls(
+            (
+                cudf.Series(point_col),
+                cudf.Series(),
+                cudf.Series(),
+                cudf.Series(),
+            ),
+            meta,
+        )

python/cuspatial/cuspatial/core/_column/geometa.py

@@ -3,11 +3,19 @@
 # This allows GeoMeta as its own init type
 from __future__ import annotations
 
+from enum import Enum
 from typing import Union
 
 import cudf
 
 
+class Feature_Enum(Enum):
+    POINT = 0
+    MULTIPOINT = 1
+    LINESTRING = 2
+    POLYGON = 3
+
+
 class GeoMeta:
     """
     Creates input_types and union_offsets for GeoColumns that are created

python/cuspatial/cuspatial/core/geoseries.py

@@ -2,7 +2,7 @@
 
 from functools import cached_property
 from numbers import Integral
-from typing import Tuple, TypeVar, Union
+from typing import Optional, Tuple, TypeVar, Union
 
 import cupy as cp
 import geopandas as gpd
@@ -22,8 +22,8 @@
 import cudf
 
 import cuspatial.io.pygeoarrow as pygeoarrow
-from cuspatial.core._column.geocolumn import GeoColumn, GeoMeta
-from cuspatial.io.geopandas_reader import Feature_Enum
+from cuspatial.core._column.geocolumn import GeoColumn
+from cuspatial.core._column.geometa import Feature_Enum, GeoMeta
 
 T = TypeVar("T", bound="GeoSeries")
 
@@ -46,14 +46,16 @@ class GeoSeries(cudf.Series):
 
     def __init__(
         self,
-        data: Union[gpd.GeoSeries, Tuple, T, pd.Series, GeoColumn],
+        data: Optional[
+            Union[gpd.GeoSeries, Tuple, T, pd.Series, GeoColumn, list]
+        ],
         index: Union[cudf.Index, pd.Index] = None,
         dtype=None,
         name=None,
         nan_as_null=True,
     ):
         # Condition data
-        if isinstance(data, pd.Series):
+        if data is None or isinstance(data, (pd.Series, list)):
             data = gpGeoSeries(data)
         # Create column
         if isinstance(data, GeoColumn):

python/cuspatial/cuspatial/core/spatial/__init__.py

@@ -26,12 +26,16 @@
     lonlat_to_cartesian,
 )
 
+from .nearest_points import pairwise_point_linestring_nearest_points
+
 __all__ = [
     "directed_hausdorff_distance",
     "haversine_distance",
     "join_quadtree_and_bounding_boxes",
     "lonlat_to_cartesian",
     "pairwise_linestring_distance",
+    "pairwise_point_linestring_distance",
+    "pairwise_point_linestring_nearest_points",
     "polygon_bounding_boxes",
     "polyline_bounding_boxes",
     "point_in_polygon",

python/cuspatial/cuspatial/core/spatial/nearest_points.py

@@ -0,0 +1,107 @@
+import cupy as cp
+
+import cuspatial._lib.nearest_points as nearest_points
+from cuspatial.core._column.geocolumn import GeoColumn
+from cuspatial.core.geodataframe import GeoDataFrame
+from cuspatial.core.geoseries import GeoSeries
+from cuspatial.utils.column_utils import (
+    contains_only_linestrings,
+    contains_only_points,
+)
+
+
+def pairwise_point_linestring_nearest_points(
+    points: GeoSeries, linestrings: GeoSeries
+) -> GeoDataFrame:
+    """Returns the nearest points between two GeoSeries of points and
+    linestrings.
+
+    Parameters
+    ----------
+    points : GeoSeries
+        A GeoSeries of points. Currently either only a series of points or
+        a series of multi-points is supported.
+    linestrings : GeoSeries
+        A GeoSeries of linestrings.
+
+    Returns
+    -------
+    GeoDataFrame
+        A GeoDataFrame with four columns.
+
+        - "point_geometry_id" contains index of the nearest point in the row.
+          If `points` consists of single-points, it is always 0.
+        - "linestring_geometry_id" contains the index of the linestring in the
+          multilinestring that contains the nearest point.
+        - "segment_id" contains the index of the segment in the linestring that
+          contains the nearest point.
+        - "geometry" contains the points of the nearest
+          point on the linestring.
+    """
+
+    if len(points) != len(linestrings):
+        raise ValueError(
+            "The inputs should have the same number of geometries"
+        )
+
+    if len(points) == 0:
+        data = {
+            "point_geometry_id": [],
+            "linestring_geometry_id": [],
+            "segment_id": [],
+            "geometry": GeoSeries([]),
+        }
+        return GeoDataFrame._from_data(data)
+
+    if not contains_only_points(points):
+        raise ValueError("`points` must contain only point geometries.")
+
+    if not contains_only_linestrings(linestrings):
+        raise ValueError(
+            "`linestrings` must contain only linestring geometries."
+        )
+
+    if len(points.points.xy) > 0 and len(points.multipoints.xy) > 0:
+        raise NotImplementedError(
+            "Mixing points and multipoint geometries in `points` is not yet "
+            "supported."
+        )
+
+    points_xy = (
+        points.points.xy
+        if len(points.points.xy) > 0
+        else points.multipoints.xy
+    )
+    points_geometry_offset = (
+        None
+        if len(points.points.xy) > 0
+        else points.multipoints.geometry_offset._column
+    )
+
+    (
+        point_geometry_id,
+        linestring_geometry_id,
+        segment_id,
+        point_on_linestring_xy,
+    ) = nearest_points.pairwise_point_linestring_nearest_points(
+        points_xy._column,
+        linestrings.lines.part_offset._column,
+        linestrings.lines.xy._column,
+        points_geometry_offset,
+        linestrings.lines.geometry_offset._column,
+    )
+
+    point_on_linestring = GeoColumn._from_points_xy(point_on_linestring_xy)
+    nearest_points_on_linestring = GeoSeries(point_on_linestring)
+
+    if not point_geometry_id:
+        point_geometry_id = cp.zeros(len(points), dtype=cp.int32)
+
+    data = {
+        "point_geometry_id": point_geometry_id,
+        "linestring_geometry_id": linestring_geometry_id,
+        "segment_id": segment_id,
+        "geometry": nearest_points_on_linestring,
+    }
+
+    return GeoDataFrame._from_data(data)