Add column API for pairwise_point_polygon_distance

cpp/CMakeLists.txt

@@ -133,6 +133,7 @@ add_library(cuspatial
     src/spatial/linestring_intersection.cu
     src/spatial/point_distance.cu
     src/spatial/point_linestring_distance.cu
+    src/spatial/point_polygon_distance.cu
     src/spatial/point_linestring_nearest_points.cu
     src/spatial/sinusoidal_projection.cu
     src/trajectory/derive_trajectories.cu

cpp/include/cuspatial/detail/utility/validation.hpp

@@ -48,3 +48,31 @@
                     "Each polygon must have at least one (1) ring");                              \
   CUSPATIAL_EXPECTS(num_poly_points >= 4 * (num_poly_ring_offsets - 1),                           \
                     "Each ring must have at least four (4) vertices");
+
+/**
+ * @brief Macro for validating the data array sizes for a multipolygon.
+ *
+ * Raises an exception if any of the following are false:
+ *  - The number of multipolygon offsets is greater than zero.
+ *  - The number of polygon offsets is greater than zero.
+ *  - The number of ring offsets is greater than zero.
+ *  - There is at least one ring offset per polygon offset.
+ *  - There are at least four vertices per ring offset.
+ *
+ * MultiPolygons follow [GeoArrow data layout][1]. Offsets arrays (polygons and rings) have one more
+ * element than the number of items in the array. The last offset is always the sum of the previous
+ * offset and the size of that element. For example the last value in the ring offsets array is the
+ * last ring offset plus the number of rings in the last polygon. See
+ * [Arrow Variable-Size Binary layout](2). Note that an empty list still has one offset: {0}.
+ *
+ * Rings are assumed to be closed (closed means the first and last vertices of
+ * each ring are equal). Therefore rings must have at least 4 vertices.
+ *
+ * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md
+ * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout
+ */
+#define CUSPATIAL_EXPECTS_VALID_MULTIPOLYGON_SIZES(                                \
+  num_poly_points, num_multipoly_offsets, num_poly_offsets, num_poly_ring_offsets) \
+  CUSPATIAL_EXPECTS(num_multipoly_offsets > 0,                                     \
+                    "Multipolygon offsets must contain at least one (1) value");   \
+  CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES(num_poly_points, num_poly_offsets, num_poly_ring_offsets);

cpp/include/cuspatial/distance/point_polygon_distance.hpp

@@ -0,0 +1,49 @@
+/*
+ * Copyright (c) 2023, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <cuspatial/column/geometry_column_view.hpp>
+
+#include <cudf/column/column_view.hpp>
+
+#include <optional>
+
+namespace cuspatial {
+
+/**
+ * @ingroup distance
+ * @brief Compute pairwise (multi)point-to-(multi)polygon Cartesian distance
+ *
+ * @param multipoints Geometry column of multipoints
+ * @param multipolygons Geometry column of multipolygons
+ * @param mr Device memory resource used to allocate the returned column.
+ * @return Column of distances between each pair of input geometries, same type as input coordinate
+ * types.
+ *
+ * @throw cuspatial::logic_error if `multipoints` and `multipolygons` has different coordinate
+ * types.
+ * @throw cuspatial::logic_error if `multipoints` is not a point column and `multipolygons` is not a
+ * polygon column.
+ * @throw cuspatial::logic_error if input column sizes mismatch.
+ */
+
+std::unique_ptr<cudf::column> pairwise_point_polygon_distance(
+  geometry_column_view const& multipoints,
+  geometry_column_view const& multipolygons,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh

@@ -40,6 +40,13 @@ namespace detail {
  * @return boolean to indicate if point is inside the polygon.
  * `false` if point is on the edge of the polygon.
  *
+ * @tparam T type of coordinate
+ * @tparam PolygonRef polygon_ref type
+ * @param test_point point to test for point in polygon
+ * @param polygon polygon to test for point in polygon
+ * @return boolean to indicate if point is inside the polygon.
+ * `false` if point is on the edge of the polygon.
+ *
  * TODO: the ultimate goal of refactoring this as independent function is to remove
  * src/utility/point_in_polygon.cuh and its usage in quadtree_point_in_polygon.cu. It isn't
  * possible today without further work to refactor quadtree_point_in_polygon into header only

cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh

@@ -23,6 +23,7 @@
 
 #include <cuspatial/cuda_utils.hpp>
 #include <cuspatial/detail/iterator.hpp>
+#include <cuspatial/detail/utility/validation.hpp>
 #include <cuspatial/experimental/geometry/segment.cuh>
 #include <cuspatial/experimental/geometry_collection/multipolygon_ref.cuh>
 #include <cuspatial/traits.hpp>
@@ -32,7 +33,6 @@
 #include <optional>
 
 namespace cuspatial {
-
 using namespace detail;
 
 template <typename GeometryIterator,
@@ -104,7 +104,10 @@ multipolygon_range<GeometryIterator, PartIterator, RingIterator, VecIterator>::m
     _point_end(point_end)
 {
   static_assert(is_vec_2d<iterator_value_type<VecIterator>>(),
-                "Coordinate range must be constructed with iterators to vec_2d.");
+                "Point iterator must be iterators to floating point vec_2d types.");
+
+  CUSPATIAL_EXPECTS_VALID_MULTIPOLYGON_SIZES(
+    num_points(), num_multipolygons() + 1, num_polygons() + 1, num_rings() + 1);
 }
 
 template <typename GeometryIterator,

cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh

@@ -245,7 +245,7 @@ auto make_multilinestring_range(GeometryColumnView const& linestrings_column)
                     "Must be Linestring geometry type.");
   auto geometry_iter       = thrust::make_counting_iterator(0);
   auto const& part_offsets = linestrings_column.offsets();
-  auto const& points_xy    = linestrings_column.child().child(1);
+  auto const& points_xy = linestrings_column.child().child(1);  // Ignores x-y offset {0, 2, 4...}
 
   auto points_it = make_vec_2d_iterator(points_xy.template begin<T>());
 
@@ -276,7 +276,7 @@ auto make_multilinestring_range(GeometryColumnView const& linestrings_column)
   auto const& geometry_offsets = linestrings_column.offsets();
   auto const& parts            = linestrings_column.child();
   auto const& part_offsets     = parts.child(0);
-  auto const& points_xy        = parts.child(1).child(1);
+  auto const& points_xy        = parts.child(1).child(1);  // Ignores x-y offset {0, 2, 4...}
 
   auto points_it = make_vec_2d_iterator(points_xy.template begin<T>());
 

cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh

@@ -18,6 +18,7 @@
 
 #include <cuspatial/cuda_utils.hpp>
 #include <cuspatial/traits.hpp>
+#include <cuspatial/types.hpp>
 
 namespace cuspatial {
 
@@ -156,7 +157,7 @@ class multipoint_range {
 };
 
 /**
- * @brief Create a multilinestring_range object of from size and start iterators
+ * @brief Create a multipoint_range object of from size and start iterators
  * @ingroup ranges
  *
  * @tparam GeometryIteratorDiffType Index type of the size of the geometry array
@@ -192,7 +193,7 @@ multipoint_range<GeometryIterator, VecIterator> make_multipoint_range(
 }
 
 /**
- * @brief Create multilinestring_range object from offset and point ranges
+ * @brief Create multipoint_range object from offset and point ranges
  *
  * @tparam IntegerRange Range to integers
  * @tparam PointRange Range to points
@@ -208,6 +209,60 @@ auto make_multipoint_range(IntegerRange geometry_offsets, PointRange points)
     geometry_offsets.begin(), geometry_offsets.end(), points.begin(), points.end());
 }
 
+/**
+ * @ingroup ranges
+ * @brief Create a range object of multipoints from cuspatial::geometry_column_view.
+ * Specialization for points column.
+ *
+ * @pre points_column must be a cuspatial::geometry_column_view
+ */
+template <collection_type_id Type,
+          typename T,
+          typename IndexType,
+          CUSPATIAL_ENABLE_IF(Type == collection_type_id::SINGLE),
+          typename GeometryColumnView>
+auto make_multipoint_range(GeometryColumnView const& points_column)
+{
+  CUSPATIAL_EXPECTS(points_column.geometry_type() == geometry_type_id::POINT,
+                    "Must be POINT geometry type.");
+  auto geometry_iter    = thrust::make_counting_iterator(0);
+  auto const& points_xy = points_column.child();  // Ignores x-y offset {0, 2, 4...}
+
+  auto points_it = make_vec_2d_iterator(points_xy.template begin<T>());
+
+  return multipoint_range(geometry_iter,
+                          thrust::next(geometry_iter, points_column.size() + 1),
+                          points_it,
+                          points_it + points_xy.size() / 2);
+}
+
+/**
+ * @ingroup ranges
+ * @brief Create a range object of multipoints from cuspatial::geometry_column_view.
+ * Specialization for multipoints column.
+ *
+ * @pre multipoints_column must be a cuspatial::geometry_column_view
+ */
+template <collection_type_id Type,
+          typename T,
+          typename IndexType,
+          CUSPATIAL_ENABLE_IF(Type == collection_type_id::MULTI),
+          typename GeometryColumnView>
+auto make_multipoint_range(GeometryColumnView const& points_column)
+{
+  CUSPATIAL_EXPECTS(points_column.geometry_type() == geometry_type_id::POINT,
+                    "Must be POINT geometry type.");
+  auto const& geometry_offsets = points_column.offsets();
+  auto const& points_xy        = points_column.child().child(1);  // Ignores x-y offset {0, 2, 4...}
+
+  auto points_it = make_vec_2d_iterator(points_xy.template begin<T>());
+
+  return multipoint_range(geometry_offsets.template begin<IndexType>(),
+                          geometry_offsets.template end<IndexType>(),
+                          points_it,
+                          points_it + points_xy.size() / 2);
+};
+
 }  // namespace cuspatial
 
 #include <cuspatial/experimental/detail/ranges/multipoint_range.cuh>

cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh

@@ -152,6 +152,74 @@ class multipolygon_range {
   CUSPATIAL_HOST_DEVICE bool is_valid_segment_id(IndexType1 segment_idx, IndexType2 ring_idx);
 };
 
+/**
+ * @ingroup ranges
+ * @brief Create a range object of multipolygon from cuspatial::geometry_column_view.
+ * Specialization for polygons column.
+ *
+ * @pre polygons_column must be a cuspatial::geometry_column_view
+ */
+template <collection_type_id Type,
+          typename T,
+          typename IndexType,
+          typename GeometryColumnView,
+          CUSPATIAL_ENABLE_IF(Type == collection_type_id::SINGLE)>
+auto make_multipolygon_range(GeometryColumnView const& polygons_column)
+{
+  CUSPATIAL_EXPECTS(polygons_column.geometry_type() == geometry_type_id::POLYGON,
+                    "Must be polygon geometry type.");
+  auto geometry_iter       = thrust::make_counting_iterator(0);
+  auto const& part_offsets = polygons_column.offsets();
+  auto const& ring_offsets = polygons_column.child().child(0);
+  auto const& points_xy =
+    polygons_column.child().child(1).child(1);  // Ignores x-y offset {0, 2, 4...}
+
+  auto points_it = make_vec_2d_iterator(points_xy.template begin<T>());
+
+  return multipolygon_range(geometry_iter,
+                            geometry_iter + part_offsets.size(),
+                            part_offsets.template begin<IndexType>(),
+                            part_offsets.template end<IndexType>(),
+                            ring_offsets.template begin<IndexType>(),
+                            ring_offsets.template end<IndexType>(),
+                            points_it,
+                            points_it + points_xy.size() / 2);
+}
+
+/**
+ * @ingroup ranges
+ * @brief Create a range object of multipolygon from cuspatial::geometry_column_view.
+ * Specialization for multipolygons column.
+ *
+ * @pre polygon_column must be a cuspatial::geometry_column_view
+ */
+template <collection_type_id Type,
+          typename T,
+          typename IndexType,
+          CUSPATIAL_ENABLE_IF(Type == collection_type_id::MULTI),
+          typename GeometryColumnView>
+auto make_multipolygon_range(GeometryColumnView const& polygons_column)
+{
+  CUSPATIAL_EXPECTS(polygons_column.geometry_type() == geometry_type_id::POLYGON,
+                    "Must be polygon geometry type.");
+  auto const& geometry_offsets = polygons_column.offsets();
+  auto const& part_offsets     = polygons_column.child().child(0);
+  auto const& ring_offsets     = polygons_column.child().child(1).child(0);
+  auto const& points_xy =
+    polygons_column.child().child(1).child(1).child(1);  // Ignores x-y offset {0, 2, 4...}
+
+  auto points_it = make_vec_2d_iterator(points_xy.template begin<T>());
+
+  return multipolygon_range(geometry_offsets.template begin<IndexType>(),
+                            geometry_offsets.template end<IndexType>(),
+                            part_offsets.template begin<IndexType>(),
+                            part_offsets.template end<IndexType>(),
+                            ring_offsets.template begin<IndexType>(),
+                            ring_offsets.template end<IndexType>(),
+                            points_it,
+                            points_it + points_xy.size() / 2);
+};
+
 }  // namespace cuspatial
 
 #include <cuspatial/experimental/detail/ranges/multipolygon_range.cuh>