Add C++ Column API and Python API for `pairwise_linestring_intersecti…

…on` (#862)

This PR adds column/python API for `pairwise_linestring_intersection`. Closes #648 

This PR is also the first attempt of #261 , establishing `geometry_column_view` that inherits `cudf::column` and adds additional information to identify the underlying nested levels. Note that it cannot directly inherit from `cudf::lists::list_column_view`, because a point array isn't a `List<T>` array.

Also note that the type codes defined in libcuspatial and cuspatial-python may be different, so there is a type code mapping after computing the result from libcuspaital.

Finally, since `List<Union>` is currently unsupported in libcudf, the python API returns the geometries in two parts: `offset + GeoColumn`.

Authors:
  - Michael Wang (https://github.com/isVoid)

Approvers:
  - H. Thomson Comer (https://github.com/thomcom)
  - Mark Harris (https://github.com/harrism)

URL: #862

cpp/CMakeLists.txt

@@ -1,5 +1,5 @@
 #=============================================================================
-# Copyright (c) 2019-2022, NVIDIA CORPORATION.
+# Copyright (c) 2019-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.
@@ -109,6 +109,7 @@ include(cmake/thirdparty/CUSPATIAL_GetCUDF.cmake)
 # - library targets -------------------------------------------------------------------------------
 
 add_library(cuspatial
+    src/column/geometry_column_view.cpp
     src/indexing/construction/point_quadtree.cu
     src/interpolate/cubic_spline.cu
     src/join/quadtree_point_in_polygon.cu
@@ -122,6 +123,7 @@ add_library(cuspatial
     src/spatial/haversine.cu
     src/spatial/hausdorff.cu
     src/spatial/linestring_distance.cu
+    src/spatial/linestring_intersection.cu
     src/spatial/point_distance.cu
     src/spatial/point_linestring_distance.cu
     src/spatial/point_linestring_nearest_points.cu

cpp/include/cuspatial/column/geometry_column_view.hpp

@@ -0,0 +1,65 @@
+/*
+ * 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/experimental/ranges/range.cuh>
+#include <cuspatial/types.hpp>
+
+#include <cudf/lists/lists_column_view.hpp>
+#include <cudf/types.hpp>
+
+namespace cuspatial {
+
+/**
+ * @ingroup cuspatial_types
+ * @brief A non-owning, immutable view of a geometry column.
+ *
+ * A geometry column is GeoArrow compliant, except that the data type for
+ * the coordinates is List<T>, instead of FixedSizeList<T>[n_dim]. This is
+ * because libcudf does not support FixedSizeList type. Currently, an
+ * even sequence (0, 2, 4, ...) is used for the offsets of the coordinate
+ * column.
+ */
+class geometry_column_view : private cudf::lists_column_view {
+ public:
+  geometry_column_view(cudf::column_view const& column,
+                       collection_type_id collection_type,
+                       geometry_type_id geometry_type);
+  geometry_column_view(geometry_column_view&&)      = default;
+  geometry_column_view(const geometry_column_view&) = default;
+  ~geometry_column_view()                           = default;
+
+  geometry_column_view& operator=(geometry_column_view const&) = default;
+
+  geometry_column_view& operator=(geometry_column_view&&) = default;
+
+  geometry_type_id geometry_type() const { return _geometry_type; }
+
+  collection_type_id collection_type() const { return _collection_type; }
+
+  cudf::data_type coordinate_type() const;
+
+  using cudf::lists_column_view::child;
+  using cudf::lists_column_view::offsets;
+  using cudf::lists_column_view::size;
+
+ protected:
+  collection_type_id _collection_type;
+  geometry_type_id _geometry_type;
+};
+
+}  // namespace cuspatial

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2022, NVIDIA CORPORATION.
+ * Copyright (c) 2022-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.
@@ -16,13 +16,14 @@
 
 #pragma once
 
-#include <thrust/pair.h>
-
 #include <cuspatial/cuda_utils.hpp>
 #include <cuspatial/experimental/detail/ranges/enumerate_range.cuh>
 #include <cuspatial/traits.hpp>
+#include <cuspatial/types.hpp>
 #include <cuspatial/vec_2d.hpp>
 
+#include <thrust/pair.h>
+
 namespace cuspatial {
 
 /**
@@ -226,6 +227,67 @@ auto make_multilinestring_range(IntegerRange1 geometry_offsets,
                                points.end());
 }
 
+/**
+ * @ingroup ranges
+ * @brief Create a range object of multilinestring from cuspatial::geometry_column_view.
+ * Specialization for linestrings column.
+ *
+ * @pre linestrings_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_multilinestring_range(GeometryColumnView const& linestrings_column)
+{
+  CUSPATIAL_EXPECTS(linestrings_column.geometry_type() == geometry_type_id::LINESTRING,
+                    "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 points_it = make_vec_2d_iterator(points_xy.template begin<T>());
+
+  return multilinestring_range(geometry_iter,
+                               geometry_iter + part_offsets.size(),
+                               part_offsets.template begin<IndexType>(),
+                               part_offsets.template end<IndexType>(),
+                               points_it,
+                               points_it + points_xy.size() / 2);
+}
+
+/**
+ * @ingroup ranges
+ * @brief Create a range object of multilinestring from cuspatial::geometry_column_view.
+ * Specialization for multilinestrings column.
+ *
+ * @pre linestring_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_multilinestring_range(GeometryColumnView const& linestrings_column)
+{
+  CUSPATIAL_EXPECTS(linestrings_column.geometry_type() == geometry_type_id::LINESTRING,
+                    "Must be Linestring geometry type.");
+  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 points_it = make_vec_2d_iterator(points_xy.template begin<T>());
+
+  return multilinestring_range(geometry_offsets.template begin<IndexType>(),
+                               geometry_offsets.template end<IndexType>(),
+                               part_offsets.template begin<IndexType>(),
+                               part_offsets.template end<IndexType>(),
+                               points_it,
+                               points_it + points_xy.size() / 2);
+};
+
 }  // namespace cuspatial
 
 #include <cuspatial/experimental/detail/ranges/multilinestring_range.cuh>

cpp/include/cuspatial/linestring_intersection.hpp

@@ -0,0 +1,59 @@
+/*
+ * 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.hpp>
+
+#include <rmm/mr/device/device_memory_resource.hpp>
+
+namespace cuspatial {
+/**
+ * @brief Result of linestring intersections
+ *
+ * Owning object to hold the result of linestring intersections.
+ * The results are modeled after arrow type List<Union<Point, LineString>>,
+ * with additional information about the indices where the intersection locates.
+ */
+struct linestring_intersection_column_result {
+  /// List offsets to the union column
+  std::unique_ptr<cudf::column> geometry_collection_offset;
+
+  /// Union Column Results
+  std::unique_ptr<cudf::column> types_buffer;
+  std::unique_ptr<cudf::column> offset_buffer;
+
+  /// Child 0: Point Results as List Type Column
+  std::unique_ptr<cudf::column> points;
+
+  /// Child 1: Segment Results as List Type Column
+  std::unique_ptr<cudf::column> segments;
+
+  /// Look-back Indices
+  std::unique_ptr<cudf::column> lhs_linestring_id;
+  std::unique_ptr<cudf::column> lhs_segment_id;
+  std::unique_ptr<cudf::column> rhs_linestring_id;
+  std::unique_ptr<cudf::column> rhs_segment_id;
+};
+
+linestring_intersection_column_result pairwise_linestring_intersection(
+  geometry_column_view const& multilinestrings_lhs,
+  geometry_column_view const& multilinestrings_rhs,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
+
+}  // namespace cuspatial

cpp/include/cuspatial/traits.hpp

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2022, NVIDIA CORPORATION.
+ * Copyright (c) 2022-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.
@@ -26,6 +26,19 @@
 
 namespace cuspatial {
 
+/**
+ * @brief Convenience macro for SFINAE as an unnamed template parameter.
+ *
+ * Example:
+ * \code{cpp}
+ * // This function will participate in overload resolution only if T is an integral type
+ * template <typename T, CUSPATIAL_ENABLE_IF(std::is_integral_v<T> )>
+ * void foo();
+ * \endcode
+ *
+ */
+#define CUSPATIAL_ENABLE_IF(...) typename std::enable_if_t<(__VA_ARGS__)>* = nullptr
+
 /**
  * @internal
  * @brief returns true if all types Ts... are the same as T.

cpp/include/cuspatial/types.hpp

@@ -0,0 +1,33 @@
+/*
+ * 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 <cstdint>
+
+namespace cuspatial {
+
+/**
+ * @brief The underlying geometry type of a geometry_column_view.
+ */
+enum class geometry_type_id : uint8_t { POINT, LINESTRING, POLYGON };
+
+/**
+ * @brief The underlying collection type of a geometry_column_view.
+ */
+enum class collection_type_id : uint8_t { SINGLE, MULTI };
+
+}  // namespace cuspatial

cpp/include/doxygen_groups.h

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2022, NVIDIA CORPORATION.
+ * Copyright (c) 2022-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.
@@ -101,7 +101,9 @@
  * @defgroup cuspatial_types Types
  * @{
  *      @brief Type declarations for cuspatial
+ *      @file types.hpp
  *      @file vec_2d.hpp
+ *      @file geometry_column_view.hpp
  *
  *      @defgroup type_factories Factory Methods
  *      @{

cpp/src/column/geometry_column_view.cpp

@@ -0,0 +1,54 @@
+/*
+ * 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.
+ */
+
+#include <cuspatial/column/geometry_column_view.hpp>
+#include <cuspatial/error.hpp>
+#include <cuspatial/types.hpp>
+
+#include <cudf/column/column_view.hpp>
+#include <cudf/lists/lists_column_view.hpp>
+#include <cudf/types.hpp>
+
+namespace cuspatial {
+
+namespace {
+
+/**
+ * @brief Retrieve the leaf column type of the LIST column.
+ * @internal
+ *
+ * @param column The LIST column to retrieve leaf column type.
+ * @return cudf type of the leaf column
+ */
+cudf::data_type leaf_data_type(cudf::column_view const& column)
+{
+  if (column.type() != cudf::data_type{cudf::type_id::LIST}) return column.type();
+  return leaf_data_type(column.child(cudf::lists_column_view::child_column_index));
+}
+
+}  // namespace
+
+geometry_column_view::geometry_column_view(cudf::column_view const& column,
+                                           collection_type_id collection_type,
+                                           geometry_type_id geometry_type)
+  : cudf::lists_column_view(column),
+    _collection_type(collection_type),
+    _geometry_type(geometry_type)
+{
+}
+
+cudf::data_type geometry_column_view::coordinate_type() const { return leaf_data_type(child()); }
+}  // namespace cuspatial