diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 1b6afb223..76d3e4198 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -134,6 +134,8 @@ add_library(cuspatial src/spatial/point_distance.cu src/spatial/point_linestring_distance.cu src/spatial/point_polygon_distance.cu + src/spatial/linestring_polygon_distance.cu + src/spatial/polygon_distance.cu src/spatial/point_linestring_nearest_points.cu src/spatial/sinusoidal_projection.cu src/trajectory/derive_trajectories.cu diff --git a/cpp/include/cuspatial/assert.cuh b/cpp/include/cuspatial/assert.cuh new file mode 100644 index 000000000..c5e13d6f4 --- /dev/null +++ b/cpp/include/cuspatial/assert.cuh @@ -0,0 +1,35 @@ +/* + * 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 + +/** + * @brief `assert`-like macro for device code + * + * This is effectively the same as the standard `assert` macro, except it + * relies on the `__PRETTY_FUNCTION__` macro which is specific to GCC and Clang + * to produce better assert messages. + */ +#if !defined(NDEBUG) && defined(__CUDA_ARCH__) && (defined(__clang__) || defined(__GNUC__)) +#define __ASSERT_STR_HELPER(x) #x +#define cuspatial_assert(e) \ + ((e) ? static_cast(0) \ + : __assert_fail(__ASSERT_STR_HELPER(e), __FILE__, __LINE__, __PRETTY_FUNCTION__)) +#else +#define cuspatial_assert(e) (static_cast(0)) +#endif diff --git a/cpp/include/cuspatial/detail/utility/validation.hpp b/cpp/include/cuspatial/detail/utility/validation.hpp index 7fe5742ea..1a4bd2001 100644 --- a/cpp/include/cuspatial/detail/utility/validation.hpp +++ b/cpp/include/cuspatial/detail/utility/validation.hpp @@ -35,10 +35,10 @@ * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout */ #define CUSPATIAL_EXPECTS_VALID_LINESTRING_SIZES(num_linestring_points, num_linestring_offsets) \ - CUSPATIAL_EXPECTS(num_linestring_offsets > 0, \ - "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_EXPECTS(num_linestring_points >= 2 * (num_linestring_offsets - 1), \ - "Each linestring must have at least two vertices"); + CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_offsets > 0, \ + "Polygon offsets must contain at least one (1) value"); \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_points >= 2 * (num_linestring_offsets - 1), \ + "Each linestring must have at least two vertices"); /** * @brief Macro for validating the data array sizes for multilinestrings. @@ -57,10 +57,10 @@ * [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_MULTILINESTRING_SIZES( \ - num_linestring_points, num_multilinestring_offsets, num_linestring_offsets) \ - CUSPATIAL_EXPECTS(num_multilinestring_offsets > 0, \ - "Multilinestring offsets must contain at least one (1) value"); \ +#define CUSPATIAL_EXPECTS_VALID_MULTILINESTRING_SIZES( \ + num_linestring_points, num_multilinestring_offsets, num_linestring_offsets) \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_multilinestring_offsets > 0, \ + "Multilinestring offsets must contain at least one (1) value"); \ CUSPATIAL_EXPECTS_VALID_LINESTRING_SIZES(num_linestring_points, num_linestring_offsets); /** @@ -84,15 +84,16 @@ * [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_POLYGON_SIZES( \ - num_poly_points, num_poly_offsets, num_poly_ring_offsets) \ - CUSPATIAL_EXPECTS(num_poly_offsets > 0, "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_EXPECTS(num_poly_ring_offsets > 0, \ - "Polygon ring offsets must contain at least one (1) value"); \ - CUSPATIAL_EXPECTS(num_poly_ring_offsets >= num_poly_offsets, \ - "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"); +#define CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES( \ + num_poly_points, num_poly_offsets, num_poly_ring_offsets) \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_offsets > 0, \ + "Polygon offsets must contain at least one (1) value"); \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets > 0, \ + "Polygon ring offsets must contain at least one (1) value"); \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets >= num_poly_offsets, \ + "Each polygon must have at least one (1) ring"); \ + CUSPATIAL_HOST_DEVICE_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. @@ -116,8 +117,8 @@ * [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"); \ +#define CUSPATIAL_EXPECTS_VALID_MULTIPOLYGON_SIZES( \ + num_poly_points, num_multipoly_offsets, num_poly_offsets, num_poly_ring_offsets) \ + CUSPATIAL_HOST_DEVICE_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); diff --git a/cpp/include/cuspatial/distance/linestring_polygon_distance.hpp b/cpp/include/cuspatial/distance/linestring_polygon_distance.hpp new file mode 100644 index 000000000..7fca8e09d --- /dev/null +++ b/cpp/include/cuspatial/distance/linestring_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 + +#include + +#include + +namespace cuspatial { + +/** + * @ingroup distance + * @brief Compute pairwise (multi)linestring-to-(multi)polygon Cartesian distance + * + * @param multilinestrings Geometry column of multilinestrings + * @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 `multilinestrings` and `multipolygons` has different coordinate + * types. + * @throw cuspatial::logic_error if `multilinestrings` is not a linestring column and + * `multipolygons` is not a polygon column. + * @throw cuspatial::logic_error if input column sizes mismatch. + */ + +std::unique_ptr pairwise_linestring_polygon_distance( + geometry_column_view const& multilinestrings, + geometry_column_view const& multipolygons, + rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource()); + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/distance/polygon_distance.hpp b/cpp/include/cuspatial/distance/polygon_distance.hpp new file mode 100644 index 000000000..ee0129992 --- /dev/null +++ b/cpp/include/cuspatial/distance/polygon_distance.hpp @@ -0,0 +1,40 @@ +/* + * 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 + +#include + +#include + +#include + +namespace cuspatial { + +/** + * @brief + * + * @param lhs + * @param rhs + * @param mr + * @return std::unique_ptr + */ +std::unique_ptr pairwise_polygon_distance( + geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource()); + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/error.hpp b/cpp/include/cuspatial/error.hpp index b1190c8ac..d2973349f 100644 --- a/cpp/include/cuspatial/error.hpp +++ b/cpp/include/cuspatial/error.hpp @@ -1,5 +1,5 @@ /* - * Copyright (c) 2020-2022, NVIDIA CORPORATION. + * Copyright (c) 2020-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,6 +16,8 @@ #pragma once +#include + #include #include @@ -76,6 +78,31 @@ struct cuda_error : public std::runtime_error { : throw cuspatial::logic_error("cuSpatial failure at: " __FILE__ \ ":" CUSPATIAL_STRINGIFY(__LINE__) ": " reason) +/**---------------------------------------------------------------------------* + * @brief Macro for checking (pre-)conditions that throws an exception when + * a condition is violated. + * + * Example usage: + * + * @code + * CUSPATIAL_HOST_DEVICE_EXPECTS(lhs->dtype == rhs->dtype, "Column type mismatch"); + * @endcode + * + * @param[in] cond Expression that evaluates to true or false + * @param[in] reason String literal description of the reason that cond is + * expected to be true + * + * (if on host) + * @throw cuspatial::logic_error if the condition evaluates to false. + * (if on device) + * program terminates and assertion error message is printed to stderr. + *---------------------------------------------------------------------------**/ +#ifndef __CUDA_ARCH__ +#define CUSPATIAL_HOST_DEVICE_EXPECTS(cond, reason) CUSPATIAL_EXPECTS(cond, reason) +#else +#define CUSPATIAL_HOST_DEVICE_EXPECTS(cond, reason) cuspatial_assert(cond&& reason) +#endif + /**---------------------------------------------------------------------------* * @brief Indicates that an erroneous code path has been taken. * diff --git a/cpp/include/cuspatial/experimental/detail/algorithm/linestring_distance.cuh b/cpp/include/cuspatial/experimental/detail/algorithm/linestring_distance.cuh new file mode 100644 index 000000000..31b984551 --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/algorithm/linestring_distance.cuh @@ -0,0 +1,74 @@ +/* + * 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 +#include + +#include + +#include + +namespace cuspatial { +namespace detail { + +/** + * @internal + * @brief The kernel to compute linestring to linestring distance + * + * Each thread of the kernel computes the distance between a segment in a linestring in pair 1 to a + * linestring in pair 2. For a segment in pair 1, the linestring index is looked up from the offset + * array and mapped to the linestring in the pair 2. The segment is then computed with all segments + * in the corresponding linestring in pair 2. This forms a local minima of the shortest distance, + * which is then combined with other segment results via an atomic operation to form the global + * minimum distance between the linestrings. + * + * `intersects` is an optional pointer to a boolean range where the `i`th element indicates the + * `i`th output should be set to 0 and bypass distance computation. This argument is optional, if + * set to nullopt, no distance computation will be bypassed. + */ +template +__global__ void linestring_distance(MultiLinestringRange1 multilinestrings1, + MultiLinestringRange2 multilinestrings2, + thrust::optional intersects, + OutputIt distances_first) +{ + using T = typename MultiLinestringRange1::element_t; + + for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multilinestrings1.num_points(); + idx += gridDim.x * blockDim.x) { + auto const part_idx = multilinestrings1.part_idx_from_point_idx(idx); + if (!multilinestrings1.is_valid_segment_id(idx, part_idx)) continue; + auto const geometry_idx = multilinestrings1.geometry_idx_from_part_idx(part_idx); + + if (intersects.has_value() && intersects.value()[geometry_idx]) { + distances_first[geometry_idx] = 0; + continue; + } + + auto [a, b] = multilinestrings1.segment(idx); + T min_distance_squared = std::numeric_limits::max(); + + for (auto const& linestring2 : multilinestrings2[geometry_idx]) { + for (auto [c, d] : linestring2) { + min_distance_squared = min(min_distance_squared, squared_segment_distance(a, b, c, d)); + } + } + atomicMin(&distances_first[geometry_idx], static_cast(sqrt(min_distance_squared))); + } +} + +} // namespace detail +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/distance_utils.cuh b/cpp/include/cuspatial/experimental/detail/distance_utils.cuh index db07dcb7d..bf17985ac 100644 --- a/cpp/include/cuspatial/experimental/detail/distance_utils.cuh +++ b/cpp/include/cuspatial/experimental/detail/distance_utils.cuh @@ -15,6 +15,7 @@ */ #include +#include #include #include #include diff --git a/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh b/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh index 97648d74b..bb146406a 100644 --- a/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/linestring_distance.cuh @@ -16,9 +16,8 @@ #pragma once -#include -#include #include +#include #include #include @@ -26,49 +25,12 @@ #include #include +#include #include #include namespace cuspatial { -namespace detail { - -/** - * @internal - * @brief The kernel to compute linestring to linestring distance - * - * Each thread of the kernel computes the distance between a segment in a linestring in pair 1 to a - * linestring in pair 2. For a segment in pair 1, the linestring index is looked up from the offset - * array and mapped to the linestring in the pair 2. The segment is then computed with all segments - * in the corresponding linestring in pair 2. This forms a local minima of the shortest distance, - * which is then combined with other segment results via an atomic operation to form the globally - * minimum distance between the linestrings. - */ -template -__global__ void pairwise_linestring_distance_kernel(MultiLinestringRange1 multilinestrings1, - MultiLinestringRange2 multilinestrings2, - OutputIt distances_first) -{ - using T = typename MultiLinestringRange1::element_t; - - for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multilinestrings1.num_points(); - idx += gridDim.x * blockDim.x) { - auto const part_idx = multilinestrings1.part_idx_from_point_idx(idx); - if (!multilinestrings1.is_valid_segment_id(idx, part_idx)) continue; - auto const geometry_idx = multilinestrings1.geometry_idx_from_part_idx(part_idx); - auto [a, b] = multilinestrings1.segment(idx); - T min_distance_squared = std::numeric_limits::max(); - - for (auto const& linestring2 : multilinestrings2[geometry_idx]) { - for (auto [c, d] : linestring2) { - min_distance_squared = min(min_distance_squared, squared_segment_distance(a, b, c, d)); - } - } - atomicMin(&distances_first[geometry_idx], static_cast(sqrt(min_distance_squared))); - } -} - -} // namespace detail template OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, @@ -98,8 +60,8 @@ OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, std::size_t const num_blocks = (multilinestrings1.num_points() + threads_per_block - 1) / threads_per_block; - detail::pairwise_linestring_distance_kernel<<>>( - multilinestrings1, multilinestrings2, distances_first); + detail::linestring_distance<<>>( + multilinestrings1, multilinestrings2, thrust::nullopt, distances_first); CUSPATIAL_CUDA_TRY(cudaGetLastError()); return distances_first + multilinestrings1.size(); diff --git a/cpp/include/cuspatial/experimental/detail/linestring_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/linestring_polygon_distance.cuh index 8e3eaac1c..0a7a5fe13 100644 --- a/cpp/include/cuspatial/experimental/detail/linestring_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/linestring_polygon_distance.cuh @@ -16,6 +16,8 @@ #pragma once +#include + #include "distance_utils.cuh" #include @@ -83,6 +85,13 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring auto geometry_id = thrust::distance(thread_bounds.begin(), it); auto local_idx = idx - *it; + // if (geometry_id == 0) + // printf("idx: %d, geometry_id: %d, intersects?: %d, local_idx: %d\n", + // static_cast(idx), + // static_cast(geometry_id), + // static_cast(intersects[geometry_id]), + // static_cast(local_idx)); + if (intersects[geometry_id]) { distances[geometry_id] = 0.0f; continue; @@ -98,9 +107,29 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring auto multipolygon_segment_id = local_idx / num_segment_this_multilinestring + multipolygons_segment_offsets[geometry_id]; + // if (geometry_id == 0) + // printf( + // "multilinestring_segment_id: %d, " + // "multipolygon_segment_id: %d\n", + // static_cast(multilinestring_segment_id), + // static_cast(multipolygon_segment_id)); + auto [a, b] = multilinestrings.segment_begin()[multilinestring_segment_id]; auto [c, d] = multipolygons.segment_begin()[multipolygon_segment_id]; + auto distance = sqrt(squared_segment_distance(a, b, c, d)); + // if (geometry_id == 0) + // printf("ab: (%f, %f) -> (%f, %f), cd: (%f, %f) -> (%f, %f), dist: %f\n", + // static_cast(a.x), + // static_cast(a.y), + // static_cast(b.x), + // static_cast(b.y), + // static_cast(c.x), + // static_cast(c.y), + // static_cast(d.x), + // static_cast(d.y), + // static_cast(distance)); + atomicMin(&distances[geometry_id], sqrt(squared_segment_distance(a, b, c, d))); } }; diff --git a/cpp/include/cuspatial/experimental/detail/polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/polygon_distance.cuh new file mode 100644 index 000000000..d4053e98d --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/polygon_distance.cuh @@ -0,0 +1,93 @@ +/* + * 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 "distance_utils.cuh" +#include "linestring_distance.cuh" + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace cuspatial { + +/** + * @brief Implementation of pairwise distance between two multipolygon ranges. + * + * All points in lhs and rhs are tested for intersection its corresponding pair, + * and if any intersection is found, the distance between the two polygons is 0. + * Otherwise, the distance is the minimum distance between any two segments in the + * multipolygon pair. + */ +template +OutputIt pairwise_polygon_distance(MultipolygonRangeA lhs, + MultipolygonRangeB rhs, + OutputIt distances_first, + rmm::cuda_stream_view stream) +{ + using T = typename MultipolygonRangeA::element_t; + + CUSPATIAL_EXPECTS(lhs.size() == rhs.size(), "Must have the same number of input rows."); + + if (lhs.size() == 0) return distances_first; + + auto lhs_as_multipoints = lhs.as_multipoint_range(); + auto rhs_as_multipoints = rhs.as_multipoint_range(); + + auto intersects = [&]() { + auto lhs_in_rhs = point_polygon_intersects(lhs_as_multipoints, rhs, stream); + auto rhs_in_lhs = point_polygon_intersects(rhs_as_multipoints, lhs, stream); + + rmm::device_uvector intersects(lhs_in_rhs.size(), stream); + thrust::transform(rmm::exec_policy(stream), + lhs_in_rhs.begin(), + lhs_in_rhs.end(), + rhs_in_lhs.begin(), + intersects.begin(), + thrust::logical_or{}); + return intersects; + }(); + + auto lhs_as_multilinestrings = lhs.as_multilinestring_range(); + auto rhs_as_multilinestrings = rhs.as_multilinestring_range(); + + thrust::fill(rmm::exec_policy(stream), + distances_first, + distances_first + lhs.size(), + std::numeric_limits::max()); + + auto [threads_per_block, num_blocks] = grid_1d(lhs.num_points()); + + detail::linestring_distance<<>>( + lhs_as_multilinestrings, rhs_as_multilinestrings, intersects.begin(), distances_first); + + CUSPATIAL_CUDA_TRY(cudaGetLastError()); + return distances_first + lhs.size(); +} + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh index 13814708e..5b1ab4ec8 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh @@ -54,10 +54,11 @@ struct to_multipoint_functor { } // namespace detail template -multipoint_range::multipoint_range(GeometryIterator geometry_begin, - GeometryIterator geometry_end, - VecIterator points_begin, - VecIterator points_end) +CUSPATIAL_HOST_DEVICE multipoint_range::multipoint_range( + GeometryIterator geometry_begin, + GeometryIterator geometry_end, + VecIterator points_begin, + VecIterator points_end) : _geometry_begin(geometry_begin), _geometry_end(geometry_end), _points_begin(points_begin), diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh index ec6e466ce..2c563f915 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -441,4 +442,39 @@ multipolygon_range:: return point_idx == _ring_begin[_part_begin[_geometry_begin[geometry_idx]]]; } +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::as_multipoint_range() +{ + auto multipoint_geometry_it = thrust::make_permutation_iterator( + _ring_begin, thrust::make_permutation_iterator(_part_begin, _geometry_begin)); + + return multipoint_range{multipoint_geometry_it, + multipoint_geometry_it + thrust::distance(_geometry_begin, _geometry_end), + _point_begin, + _point_end}; +} + +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range:: + as_multilinestring_range() +{ + auto multilinestring_geometry_it = + thrust::make_permutation_iterator(_part_begin, _geometry_begin); + return multilinestring_range{ + multilinestring_geometry_it, + multilinestring_geometry_it + thrust::distance(_geometry_begin, _geometry_end), + _ring_begin, + _ring_end, + _point_begin, + _point_end}; +} + } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/polygon_distance.cuh b/cpp/include/cuspatial/experimental/polygon_distance.cuh new file mode 100644 index 000000000..63191aa71 --- /dev/null +++ b/cpp/include/cuspatial/experimental/polygon_distance.cuh @@ -0,0 +1,47 @@ +/* + * 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 + +namespace cuspatial { + +/** + * @ingroup distance + * @brief Computes pairwise multipolygon to multipolygon distance + * + * @tparam MultiPolygonRangeA An instance of template type `multipolygon_range` + * @tparam MultiPolygonRangeB An instance of template type `multipolygon_range` + * @tparam OutputIt iterator type for output array. Must meet the requirements of [LRAI](LinkLRAI). + * Must be an iterator to type convertible from floating points. + * + * @param lhs The first multipolygon range to compute distance from + * @param rhs The second multipolygon range to compute distance to + * @param stream The CUDA stream on which to perform computations + * @return Output Iterator past the last distance computed + * + * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator + * "LegacyRandomAccessIterator" + */ +template +OutputIt pairwise_polygon_distance(MultipolygonRangeA lhs, + MultipolygonRangeB rhs, + OutputIt distances_first, + rmm::cuda_stream_view stream = rmm::cuda_stream_default); +} // namespace cuspatial + +#include diff --git a/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh b/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh index 3637c327e..3caa568e2 100644 --- a/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multilinestring_range.cuh @@ -59,12 +59,12 @@ class multilinestring_range { using point_t = iterator_value_type; using element_t = iterator_vec_base_type; - multilinestring_range(GeometryIterator geometry_begin, - GeometryIterator geometry_end, - PartIterator part_begin, - PartIterator part_end, - VecIterator points_begin, - VecIterator points_end); + CUSPATIAL_HOST_DEVICE multilinestring_range(GeometryIterator geometry_begin, + GeometryIterator geometry_end, + PartIterator part_begin, + PartIterator part_end, + VecIterator points_begin, + VecIterator points_end); /// Return the number of multilinestrings in the array. CUSPATIAL_HOST_DEVICE auto size() { return num_multilinestrings(); } diff --git a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh index 7117d28cc..82323e9ed 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh @@ -57,10 +57,10 @@ class multipoint_range { /** * @brief Construct a new multipoint array object */ - multipoint_range(GeometryIterator geometry_begin, - GeometryIterator geometry_end, - VecIterator points_begin, - VecIterator points_end); + CUSPATIAL_HOST_DEVICE multipoint_range(GeometryIterator geometry_begin, + GeometryIterator geometry_end, + VecIterator points_begin, + VecIterator points_end); /** * @brief Returns the number of multipoints in the array. */ diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh index 96129145f..6c7ef781d 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -62,7 +62,9 @@ class multipolygon_range { using ring_it_t = RingIterator; using point_it_t = VecIterator; using point_t = iterator_value_type; - using element_t = iterator_vec_base_type; + + using index_t = iterator_value_type; + using element_t = iterator_vec_base_type; int64_t static constexpr INVALID_INDEX = -1; @@ -180,6 +182,15 @@ class multipolygon_range { /// Returns an iterator to the end of the segment CUSPATIAL_HOST_DEVICE auto segment_end(); + /// Range Casting + + /// Cast the range of multipolygons as a range of multipoints, ignoring all edge connections and + /// ring relationships. + CUSPATIAL_HOST_DEVICE auto as_multipoint_range(); + + /// Cast the range of multipolygons as a range of multilinestrings, ignoring ring relationships. + CUSPATIAL_HOST_DEVICE auto as_multilinestring_range(); + protected: GeometryIterator _geometry_begin; GeometryIterator _geometry_end; diff --git a/cpp/include/cuspatial_test/base_fixture.hpp b/cpp/include/cuspatial_test/base_fixture.hpp index ae92c9fd3..941bee951 100644 --- a/cpp/include/cuspatial_test/base_fixture.hpp +++ b/cpp/include/cuspatial_test/base_fixture.hpp @@ -67,7 +67,9 @@ class BaseFixture : public RMMResourceMixin, public ::testing::Test { * class MyTest : public cuspatial::test::BaseFixtureWithParam {}; * * TEST_P(MyTest, TestParamterGet) { - * auto [a, b, c] = GetParam(); + * auto a = std::get<0>(GetParam()); + * auto b = std::get<1>(GetParam()); + * auto c = std::get<2>(GetParam()); * ... * } * diff --git a/cpp/include/cuspatial_test/test_util.cuh b/cpp/include/cuspatial_test/test_util.cuh index e30089ba4..2312fb85f 100644 --- a/cpp/include/cuspatial_test/test_util.cuh +++ b/cpp/include/cuspatial_test/test_util.cuh @@ -100,7 +100,14 @@ void print_device_range(Iter begin, } /** - * @brief + * @brief Print a device vector. + * + * @note Copies the device vector to host before printing. + * + * @tparam Vector The device vector type + * @param vec The device vector to print + * @param pre String to print before the device vector + * @param post String to print after the device vector */ template void print_device_vector(Vector const& vec, std::string_view pre = "", std::string_view post = "\n") diff --git a/cpp/include/cuspatial_test/vector_equality.hpp b/cpp/include/cuspatial_test/vector_equality.hpp index 204a1589d..b4ce1ca68 100644 --- a/cpp/include/cuspatial_test/vector_equality.hpp +++ b/cpp/include/cuspatial_test/vector_equality.hpp @@ -238,6 +238,39 @@ void expect_vec_2d_pair_equivalent(PairVector1 const& expected, PairVector2 cons cuspatial::test::expect_vec_2d_pair_equivalent(lhs, rhs); \ } while (0) +template +void expect_multilinestring_array_equivalent(Array1& lhs, Array2& rhs) +{ + auto [lhs_geometry_offset, lhs_part_offset, lhs_coordinates] = lhs.release(); + auto [rhs_geometry_offset, rhs_part_offset, rhs_coordinates] = rhs.release(); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(lhs_geometry_offset, rhs_geometry_offset); + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(lhs_part_offset, rhs_part_offset); + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(lhs_coordinates, rhs_coordinates); +} + +#define CUSPATIAL_EXPECT_MULTILINESTRING_ARRAY_EQUIVALENT(lhs, rhs) \ + do { \ + SCOPED_TRACE(" <-- line of failure\n"); \ + cuspatial::test::expect_multilinestring_array_equivalent(lhs, rhs); \ + } while (0) + +template +void expect_multipoint_array_equivalent(Array1& lhs, Array2& rhs) +{ + auto [lhs_geometry_offset, lhs_coordinates] = lhs.release(); + auto [rhs_geometry_offset, rhs_coordinates] = rhs.release(); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(lhs_geometry_offset, rhs_geometry_offset); + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(lhs_coordinates, rhs_coordinates); +} + +#define CUSPATIAL_EXPECT_MULTIPOINT_ARRAY_EQUIVALENT(lhs, rhs) \ + do { \ + SCOPED_TRACE(" <-- line of failure\n"); \ + cuspatial::test::expect_multipoint_array_equivalent(lhs, rhs); \ + } while (0) + #define CUSPATIAL_RUN_TEST(FUNC, ...) \ do { \ SCOPED_TRACE(" <-- line of failure\n"); \ diff --git a/cpp/include/cuspatial_test/vector_factories.cuh b/cpp/include/cuspatial_test/vector_factories.cuh index 0e5ae930a..892f73a7a 100644 --- a/cpp/include/cuspatial_test/vector_factories.cuh +++ b/cpp/include/cuspatial_test/vector_factories.cuh @@ -240,12 +240,43 @@ class multilinestring_array { _coordinate_array.end()); } + auto release() + { + return std::tuple{std::move(_geometry_offset_array), + std::move(_part_offset_array), + std::move(_coordinate_array)}; + } + protected: GeometryArray _geometry_offset_array; PartArray _part_offset_array; CoordinateArray _coordinate_array; }; +/** + * @brief Construct an owning object of a multilinestring array from ranges + * + * @param geometry_inl Range of geometry offsets + * @param part_inl Range of part offsets + * @param coord_inl Ramge of coordinate + * @return multilinestring array object + */ +template +auto make_multilinestring_array(IndexRangeA geometry_inl, + IndexRangeB part_inl, + CoordRange coord_inl) +{ + using CoordType = typename CoordRange::value_type; + using DeviceIndexVector = thrust::device_vector; + using DeviceCoordVector = thrust::device_vector; + + return multilinestring_array( + make_device_vector(geometry_inl), make_device_vector(part_inl), make_device_vector(coord_inl)); +} + /** * @brief Construct an owning object of a multilinestring array from initializer lists * @@ -299,6 +330,16 @@ class multipoint_array { CoordinateArray _coordinates; }; +/** + * @brief Factory method to construct multipoint array from ranges of geometry offsets and + * coordinates + */ +template +auto make_multipoints_array(GeometryRange geometry_inl, CoordRange coordinates_inl) +{ + return multipoint_array{make_device_vector(geometry_inl), make_device_vector(coordinates_inl)}; +} + /** * @brief Factory method to construct multipoint array from initializer list of multipoints. * diff --git a/cpp/src/spatial/linestring_polygon_distance.cu b/cpp/src/spatial/linestring_polygon_distance.cu new file mode 100644 index 000000000..ada232ea1 --- /dev/null +++ b/cpp/src/spatial/linestring_polygon_distance.cu @@ -0,0 +1,127 @@ +/* + * 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 "../utility/iterator.hpp" +#include "../utility/multi_geometry_dispatch.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace cuspatial { + +namespace detail { + +namespace { + +template +struct pairwise_linestring_polygon_distance_impl { + using SizeType = cudf::device_span::size_type; + + template )> + std::unique_ptr operator()(geometry_column_view const& multilinestrings, + geometry_column_view const& multipolygons, + rmm::cuda_stream_view stream, + rmm::mr::device_memory_resource* mr) + { + auto multilinestrings_range = + make_multilinestring_range(multilinestrings); + auto multipolygons_range = + make_multipolygon_range(multipolygons); + + auto output = cudf::make_numeric_column(multilinestrings.coordinate_type(), + multilinestrings.size(), + cudf::mask_state::UNALLOCATED, + stream, + mr); + + cuspatial::pairwise_linestring_polygon_distance( + multilinestrings_range, multipolygons_range, output->mutable_view().begin(), stream); + return output; + } + + template ), typename... Args> + std::unique_ptr operator()(Args&&...) + + { + CUSPATIAL_FAIL("linestring-polygon distance API only supports floating point coordinates."); + } +}; + +} // namespace + +template +struct pairwise_linestring_polygon_distance { + std::unique_ptr operator()(geometry_column_view const& multilinestrings, + geometry_column_view const& multipolygons, + rmm::cuda_stream_view stream, + rmm::mr::device_memory_resource* mr) + { + return cudf::type_dispatcher( + multilinestrings.coordinate_type(), + pairwise_linestring_polygon_distance_impl{}, + multilinestrings, + multipolygons, + stream, + mr); + } +}; + +} // namespace detail + +std::unique_ptr pairwise_linestring_polygon_distance( + geometry_column_view const& multilinestrings, + geometry_column_view const& multipolygons, + rmm::mr::device_memory_resource* mr) +{ + CUSPATIAL_EXPECTS(multilinestrings.geometry_type() == geometry_type_id::LINESTRING && + multipolygons.geometry_type() == geometry_type_id::POLYGON, + "Unexpected input geometry types."); + + CUSPATIAL_EXPECTS(multilinestrings.coordinate_type() == multipolygons.coordinate_type(), + "Input geometries must have the same coordinate data types."); + + return multi_geometry_double_dispatch( + multilinestrings.collection_type(), + multipolygons.collection_type(), + multilinestrings, + multipolygons, + rmm::cuda_stream_default, + mr); +} + +} // namespace cuspatial diff --git a/cpp/src/spatial/polygon_distance.cu b/cpp/src/spatial/polygon_distance.cu new file mode 100644 index 000000000..65143e137 --- /dev/null +++ b/cpp/src/spatial/polygon_distance.cu @@ -0,0 +1,113 @@ +/* + * 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 "../utility/iterator.hpp" +#include "../utility/multi_geometry_dispatch.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace cuspatial { + +namespace detail { + +namespace { + +template +struct pairwise_polygon_distance_impl { + using SizeType = cudf::device_span::size_type; + + template )> + std::unique_ptr operator()(geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::cuda_stream_view stream, + rmm::mr::device_memory_resource* mr) + { + auto lhs_range = make_multipolygon_range(lhs); + auto rhs_range = make_multipolygon_range(rhs); + + auto output = cudf::make_numeric_column( + lhs.coordinate_type(), lhs.size(), cudf::mask_state::UNALLOCATED, stream, mr); + + cuspatial::pairwise_polygon_distance( + lhs_range, rhs_range, output->mutable_view().begin(), stream); + return output; + } + + template ), typename... Args> + std::unique_ptr operator()(Args&&...) + + { + CUSPATIAL_FAIL("polygon distance API only supports floating point coordinates."); + } +}; + +} // namespace + +template +struct pairwise_polygon_distance { + std::unique_ptr operator()(geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::cuda_stream_view stream, + rmm::mr::device_memory_resource* mr) + { + return cudf::type_dispatcher( + lhs.coordinate_type(), + pairwise_polygon_distance_impl{}, + lhs, + rhs, + stream, + mr); + } +}; + +} // namespace detail + +std::unique_ptr pairwise_polygon_distance(geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::mr::device_memory_resource* mr) +{ + CUSPATIAL_EXPECTS(lhs.geometry_type() == geometry_type_id::POLYGON && + rhs.geometry_type() == geometry_type_id::POLYGON, + "Unexpected input geometry types."); + + CUSPATIAL_EXPECTS(lhs.coordinate_type() == rhs.coordinate_type(), + "Input geometries must have the same coordinate data types."); + + return multi_geometry_double_dispatch( + lhs.collection_type(), rhs.collection_type(), lhs, rhs, rmm::cuda_stream_default, mr); +} + +} // namespace cuspatial diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index aa37e25cb..ffa8903de 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -90,6 +90,12 @@ ConfigureTest(LINESTRING_DISTANCE_TEST ConfigureTest(POINT_POLYGON_DISTANCE_TEST spatial/point_polygon_distance_test.cpp) +ConfigureTest(LINESTRING_POLYGON_DISTANCE_TEST + spatial/linestring_polygon_distance_test.cpp) + +ConfigureTest(POLYGON_DISTANCE_TEST + spatial/polygon_distance_test.cpp) + ConfigureTest(LINESTRING_INTERSECTION_TEST spatial/linestring_intersection_test.cpp) @@ -133,6 +139,9 @@ ConfigureTest(POINT_POLYGON_DISTANCE_TEST_EXP ConfigureTest(LINESTRING_POLYGON_DISTANCE_TEST_EXP experimental/spatial/linestring_polygon_distance_test.cu) +ConfigureTest(POLYGON_DISTANCE_TEST_EXP + experimental/spatial/polygon_distance_test.cu) + ConfigureTest(HAUSDORFF_TEST_EXP experimental/spatial/hausdorff_test.cu) diff --git a/cpp/tests/experimental/range/multipolygon_range_test.cu b/cpp/tests/experimental/range/multipolygon_range_test.cu index b6b5bc340..094abd293 100644 --- a/cpp/tests/experimental/range/multipolygon_range_test.cu +++ b/cpp/tests/experimental/range/multipolygon_range_test.cu @@ -99,6 +99,56 @@ struct MultipolygonRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, d_expected); } + + void test_multipolygon_as_multilinestring( + std::initializer_list multipolygon_geometry_offset, + std::initializer_list multipolygon_part_offset, + std::initializer_list ring_offset, + std::initializer_list> multipolygon_coordinates, + std::initializer_list multilinestring_geometry_offset, + std::initializer_list multilinestring_part_offset, + std::initializer_list> multilinestring_coordinates) + { + auto multipolygon_array = make_multipolygon_array(multipolygon_geometry_offset, + multipolygon_part_offset, + ring_offset, + multipolygon_coordinates); + auto rng = multipolygon_array.range().as_multilinestring_range(); + + auto got = + make_multilinestring_array(range(rng.geometry_offsets_begin(), rng.geometry_offsets_end()), + range(rng.part_offsets_begin(), rng.part_offsets_end()), + range(rng.point_begin(), rng.point_end())); + + auto expected = make_multilinestring_array( + multilinestring_geometry_offset, multilinestring_part_offset, multilinestring_coordinates); + + CUSPATIAL_EXPECT_MULTILINESTRING_ARRAY_EQUIVALENT(expected, got); + } + + void test_multipolygon_as_multipoint( + std::initializer_list multipolygon_geometry_offset, + std::initializer_list multipolygon_part_offset, + std::initializer_list ring_offset, + std::initializer_list> multipolygon_coordinates, + std::initializer_list multipoint_geometry_offset, + std::initializer_list> multipoint_coordinates) + { + auto multipolygon_array = make_multipolygon_array(multipolygon_geometry_offset, + multipolygon_part_offset, + ring_offset, + multipolygon_coordinates); + auto rng = multipolygon_array.range().as_multipoint_range(); + + auto got = make_multipoints_array(range(rng.offsets_begin(), rng.offsets_end()), + range(rng.point_begin(), rng.point_end())); + + auto expected = make_multipoints_array( + range(multipoint_geometry_offset.begin(), multipoint_geometry_offset.end()), + range(multipoint_coordinates.begin(), multipoint_coordinates.end())); + + CUSPATIAL_EXPECT_MULTIPOINT_ARRAY_EQUIVALENT(expected, got); + } }; TYPED_TEST_CASE(MultipolygonRangeTest, FloatingPointTypes); @@ -338,3 +388,160 @@ TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmp {0, 0}}, {6, 3}); } + +TYPED_TEST(MultipolygonRangeTest, MultipolygonAsMultilinestring1) +{ + CUSPATIAL_RUN_TEST(this->test_multipolygon_as_multilinestring, + {0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {0, 1, 2}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}); +} + +TYPED_TEST(MultipolygonRangeTest, MultipolygonAsMultilinestring2) +{ + CUSPATIAL_RUN_TEST(this->test_multipolygon_as_multilinestring, + {0, 1, 2}, + {0, 1, 3}, + {0, 4, 8, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}, + {0, 1, 3}, + {0, 4, 8, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}); +} + +TYPED_TEST(MultipolygonRangeTest, MultipolygonAsMultilinestring3) +{ + CUSPATIAL_RUN_TEST(this->test_multipolygon_as_multilinestring, + {0, 1, 2}, + {0, 2, 3}, + {0, 4, 8, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}, + {0, 2, 3}, + {0, 4, 8, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}); +} + +TYPED_TEST(MultipolygonRangeTest, MultipolygonAsMultiPoint1) +{ + CUSPATIAL_RUN_TEST(this->test_multipolygon_as_multipoint, + {0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}); +} + +TYPED_TEST(MultipolygonRangeTest, MultipolygonAsMultiPoint2) +{ + CUSPATIAL_RUN_TEST(this->test_multipolygon_as_multipoint, + {0, 1, 2}, + {0, 1, 3}, + {0, 4, 8, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}, + {0, 4, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}); +} + +TYPED_TEST(MultipolygonRangeTest, MultipolygonAsMultiPoint3) +{ + CUSPATIAL_RUN_TEST(this->test_multipolygon_as_multipoint, + {0, 1, 2}, + {0, 2, 3}, + {0, 4, 8, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}, + {0, 8, 12}, + {{0, 0}, + {1, 0}, + {1, 1}, + {0, 0}, + {10, 10}, + {11, 10}, + {11, 11}, + {10, 10}, + {20, 20}, + {21, 20}, + {21, 21}, + {20, 20}}); +} diff --git a/cpp/tests/experimental/spatial/polygon_distance_test.cu b/cpp/tests/experimental/spatial/polygon_distance_test.cu new file mode 100644 index 000000000..3d3e42b6e --- /dev/null +++ b/cpp/tests/experimental/spatial/polygon_distance_test.cu @@ -0,0 +1,482 @@ +/* + * 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 +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +using namespace cuspatial; +using namespace cuspatial::test; + +template +struct PairwisePolygonDistanceTest : BaseFixture { + template + void run_test(MultipolygonRangeA lhs, + MultipolygonRangeB rhs, + rmm::device_uvector const& expected) + { + auto got = rmm::device_uvector(lhs.size(), stream()); + auto ret = pairwise_polygon_distance(lhs, rhs, got.begin(), stream()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(expected, got); + EXPECT_EQ(thrust::distance(got.begin(), ret), expected.size()); + } + + void run(std::initializer_list lhs_multipolygon_geometry_offsets, + std::initializer_list lhs_multipolygon_part_offsets, + std::initializer_list lhs_multipolygon_ring_offsets, + std::initializer_list> lhs_multipolygon_coordinates, + std::initializer_list rhs_multipolygon_geometry_offsets, + std::initializer_list rhs_multipolygon_part_offsets, + std::initializer_list rhs_multipolygon_ring_offsets, + std::initializer_list> rhs_multipolygon_coordinates, + std::initializer_list expected) + { + auto lhs = make_multipolygon_array(lhs_multipolygon_geometry_offsets, + lhs_multipolygon_part_offsets, + lhs_multipolygon_ring_offsets, + lhs_multipolygon_coordinates); + + auto rhs = make_multipolygon_array(rhs_multipolygon_geometry_offsets, + rhs_multipolygon_part_offsets, + rhs_multipolygon_ring_offsets, + rhs_multipolygon_coordinates); + + auto lhs_range = lhs.range(); + auto rhs_range = rhs.range(); + + auto d_expected = make_device_uvector(expected, stream(), mr()); + + // Euclidean distance is symmetric + run_test(lhs_range, rhs_range, d_expected); + run_test(rhs_range, lhs_range, d_expected); + } +}; + +TYPED_TEST_CASE(PairwisePolygonDistanceTest, FloatingPointTypes); + +TYPED_TEST(PairwisePolygonDistanceTest, Empty) +{ + this->run({0}, {0}, {0}, {}, {0}, {0}, {0}, {}, {}); +} + +// Test Matrix +// One Pair: +// lhs-rhs Relationship: Disjoint, Touching, Overlapping, Contained, Within +// Holes: No, Yes +// Multipolygon: No, Yes + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonDisjointNoHole) +{ + this->run({0, 1}, + {0, 1}, + {0, 4}, + {{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{-1, 0}, {-1, 1}, {-2, 0}, {-1, 0}}, + {1}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonTouchingNoHole) +{ + this->run({0, 1}, + {0, 1}, + {0, 4}, + {{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{1, 0}, {2, 0}, {2, 1}, {1, 0}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonOverlappingNoHole) +{ + this->run({0, 1}, + {0, 1}, + {0, 4}, + {{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{0.5, 0}, {2, 0}, {2, 1}, {0.5, 0}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonContainedNoHole) +{ + this->run({0, 1}, + {0, 1}, + {0, 5}, + {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}}, + {0, 1}, + {0, 1}, + {0, 5}, + {{0.25, 0.25}, {0.75, 0.25}, {0.75, 0.75}, {0.25, 0.75}, {0.25, 0.25}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonWithinNoHole) +{ + this->run({0, 1}, + {0, 1}, + {0, 5}, + {{0.25, 0.25}, {0.75, 0.25}, {0.75, 0.75}, {0.25, 0.75}, {0.25, 0.25}}, + {0, 1}, + {0, 1}, + {0, 5}, + {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonDisjointHasHole) +{ + this->run({0, 1}, + {0, 2}, + {0, 4, 8}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {1.0, 0.75}, + {1.5, 0.75}, + {1.25, 1.0}, + {1.0, 0.75}}, + {0, 1}, + {0, 2}, + {0, 4, 8}, + {{-1.0, 0.0}, + {-1.0, -1.0}, + {-2.0, 0.0}, + {-1.0, 0.0}, + {-1.25, -0.25}, + {-1.25, -0.5}, + {-1.5, -0.25}, + {-1.25, -0.25}}, + {1}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonDisjointHasHole2) +{ + this->run({0, 1}, + {0, 2}, + {0, 5, 10}, + {{0.0, 0.0}, + {10.0, 0.0}, + {10.0, 10.0}, + {0.0, 10.0}, + {0.0, 0.0}, + {2.0, 2.0}, + {2.0, 6.0}, + {6.0, 6.0}, + {6.0, 2.0}, + {2.0, 2.0}}, + {0, 1}, + {0, 1}, + {0, 5}, + {{3.0, 3.0}, {3.0, 4.0}, {4.0, 4.0}, {4.0, 3.0}, {3.0, 3.0}}, + {1}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonTouchingHasHole) +{ + this->run({0, 1}, + {0, 2}, + {0, 4, 8}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {1.0, 0.75}, + {1.5, 0.75}, + {1.25, 1.0}, + {1.0, 0.75}}, + {0, 1}, + {0, 2}, + {0, 4, 8}, + {{2.0, 0.0}, + {3.0, 0.0}, + {3.0, 1.0}, + {2.0, 0.0}, + {2.5, 0.25}, + {2.75, 0.25}, + {2.75, 0.5}, + {2.5, 0.25}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonOverlappingHasHole) +{ + this->run({0, 1}, + {0, 2}, + {0, 5, 10}, + {{0, 0}, {4, 0}, {4, 4}, {0, 4}, {0, 0}, {2, 2}, {2, 3}, {3, 3}, {3, 2}, {2, 2}}, + {0, 1}, + {0, 2}, + {0, 4, 8}, + {{2, -1}, {5, 4}, {5, -1}, {2, -1}, {3, -0.5}, {4, 0}, {4, -0.5}, {3, -0.5}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonContainedHasHole) +{ + this->run({0, 1}, + {0, 2}, + {0, 5, 10}, + {{0, 0}, {4, 0}, {4, 4}, {0, 4}, {0, 0}, {1, 3}, {1, 1}, {3, 1}, {1, 3}, {1, 1}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{1, 3}, {3, 1}, {3, 3}, {1, 3}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairSinglePolygonWithinHasHole) +{ + this->run({0, 1}, + {0, 1}, + {0, 4}, + {{1, 3}, {3, 1}, {3, 3}, {1, 3}}, + {0, 1}, + {0, 2}, + {0, 5, 9}, + {{0, 0}, {4, 0}, {4, 4}, {0, 4}, {0, 0}, {1, 1}, {3, 1}, {1, 3}, {1, 1}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairMultiPolygonDisjointNoHole) +{ + this->run({0, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {3.0, 3.0}, + {3.0, 4.0}, + {4.0, 4.0}, + {3.0, 3.0}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{-1.0, 0.0}, {-1.0, -1.0}, {-2.0, 0.0}, {-1.0, 0.0}}, + {1}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairMultiPolygonTouchingNoHole) +{ + this->run({0, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {3.0, 3.0}, + {3.0, 4.0}, + {4.0, 4.0}, + {3.0, 3.0}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{3.0, 3.0}, {2.0, 3.0}, {2.0, 2.0}, {3.0, 3.0}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairMultiPolygonOverlappingNoHole) +{ + this->run({0, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {3.0, 3.0}, + {3.0, 4.0}, + {4.0, 4.0}, + {3.0, 3.0}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{1.0, 1.0}, {3.0, 1.0}, {3.0, 3.0}, {1.0, 1.0}}, + {0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, OnePairMultiPolygonContainedNoHole) +{ + this->run({0, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {1.0, 1.0}, + {1.0, 2.0}, + {2.0, 2.0}, + {1.0, 1.0}}, + {0, 1}, + {0, 1}, + {0, 4}, + {{0.5, 0.25}, {1.5, 0.25}, {1.5, 1.25}, {0.5, 0.25}}, + {0}); +} + +// Two Pair Tests + +TYPED_TEST(PairwisePolygonDistanceTest, TwoPairSinglePolygonNoHole) +{ + this->run({0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {3.0, 3.0}, + {3.0, 4.0}, + {4.0, 4.0}, + {3.0, 3.0}}, + {0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{3.0, 0.0}, + {4.0, 0.0}, + {4.0, 1.0}, + {3.0, 0.0}, + {3.0, 3.0}, + {3.0, 2.0}, + {2.0, 2.0}, + {3.0, 3.0}}, + {1, 0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, TwoPairSinglePolygonHasHole) +{ + this->run({0, 1, 2}, + {0, 1, 3}, + {0, 4, 8, 12}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {3.0, 3.0}, + {3.0, 4.0}, + {4.0, 4.0}, + {3.0, 3.0}, + {3.25, 3.5}, + {3.5, 3.5}, + {3.5, 3.75}, + {3.25, 3.5}}, + {0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{3.0, 0.0}, + {4.0, 0.0}, + {4.0, 1.0}, + {3.0, 0.0}, + {3.0, 3.0}, + {3.0, 2.0}, + {2.0, 2.0}, + {3.0, 3.0}}, + {1, 0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, TwoPairMultiPolygonNoHole) +{ + this->run({0, 1, 3}, + {0, 1, 2, 3}, + {0, 4, 8, 12}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + {3.0, 3.0}, + {3.0, 4.0}, + {4.0, 4.0}, + {3.0, 3.0}, + {3.0, 3.0}, + {5.0, 3.0}, + {4.0, 2.0}, + {3.0, 3.0}}, + {0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{3.0, 0.0}, + {4.0, 0.0}, + {4.0, 1.0}, + {3.0, 0.0}, + {3.0, 3.0}, + {3.0, 2.0}, + {2.0, 2.0}, + {3.0, 3.0}}, + {1, 0}); +} + +TYPED_TEST(PairwisePolygonDistanceTest, TwoPairMultiPolygonHasHole) +{ + this->run({0, 1, 3}, + {0, 1, 2, 4}, + {0, 4, 8, 12, 16}, + {{0.0, 0.0}, + {2.0, 0.0}, + {2.0, 2.0}, + {0.0, 0.0}, + + {3.0, 3.0}, + {3.0, 4.0}, + {4.0, 4.0}, + {3.0, 3.0}, + + {3.0, 3.0}, + {5.0, 3.0}, + {4.0, 2.0}, + {3.0, 3.0}, + + {3.5, 2.9}, + {4.5, 2.9}, + {4, 2.4}, + {3.5, 2.9}}, + {0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{3.0, 0.0}, + {4.0, 0.0}, + {4.0, 1.0}, + {3.0, 0.0}, + {3.0, 3.0}, + {3.0, 2.0}, + {2.0, 2.0}, + {3.0, 3.0}}, + {1, 0}); +} diff --git a/cpp/tests/spatial/linestring_polygon_distance_test.cpp b/cpp/tests/spatial/linestring_polygon_distance_test.cpp new file mode 100644 index 000000000..fdef0661b --- /dev/null +++ b/cpp/tests/spatial/linestring_polygon_distance_test.cpp @@ -0,0 +1,230 @@ +/* + * 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 +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include +#include + +using namespace cuspatial; +using namespace cuspatial::test; + +using namespace cudf; +using namespace cudf::test; + +template +struct PairwiseLinestringPolygonDistanceTest : ::testing::Test { + rmm::cuda_stream_view stream() { return cudf::get_default_stream(); } + void SetUp() + { + collection_type_id _; + std::tie(_, empty_linestring_column) = make_linestring_column({0}, {}, stream()); + std::tie(_, empty_multilinestring_column) = make_linestring_column({0}, {0}, {}, stream()); + std::tie(_, empty_polygon_column) = make_polygon_column({0}, {0}, {}, stream()); + std::tie(_, empty_multipolygon_column) = make_polygon_column({0}, {0}, {0}, {}, stream()); + } + + geometry_column_view empty_linestring() + { + return geometry_column_view( + empty_linestring_column->view(), collection_type_id::SINGLE, geometry_type_id::LINESTRING); + } + + geometry_column_view empty_multilinestring() + { + return geometry_column_view(empty_multilinestring_column->view(), + collection_type_id::MULTI, + geometry_type_id::LINESTRING); + } + + geometry_column_view empty_polygon() + { + return geometry_column_view( + empty_polygon_column->view(), collection_type_id::SINGLE, geometry_type_id::POLYGON); + } + + geometry_column_view empty_multipolygon() + { + return geometry_column_view( + empty_multipolygon_column->view(), collection_type_id::MULTI, geometry_type_id::POLYGON); + } + + void run_single(geometry_column_view linestrings, + geometry_column_view polygons, + std::initializer_list expected) + { + auto got = pairwise_linestring_polygon_distance(linestrings, polygons); + CUDF_TEST_EXPECT_COLUMNS_EQUIVALENT(*got, fixed_width_column_wrapper(expected)); + } + + std::unique_ptr empty_linestring_column; + std::unique_ptr empty_multilinestring_column; + std::unique_ptr empty_polygon_column; + std::unique_ptr empty_multipolygon_column; +}; + +struct PairwiseLinestringPolygonDistanceTestUntyped : testing::Test { + rmm::cuda_stream_view stream() { return cudf::get_default_stream(); } +}; + +using TestTypes = ::testing::Types; + +TYPED_TEST_CASE(PairwiseLinestringPolygonDistanceTest, TestTypes); + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, SingleToSingleEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_linestring(), this->empty_polygon(), {}); +}; + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, SingleToMultiEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_linestring(), this->empty_multipolygon(), {}); +}; + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, MultiToSingleEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_multilinestring(), this->empty_polygon(), {}); +}; + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, MultiToMultiEmpty) +{ + CUSPATIAL_RUN_TEST( + this->run_single, this->empty_multilinestring(), this->empty_multipolygon(), {}); +}; + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, SingleToSingleOnePair) +{ + using T = TypeParam; + + auto [ptype, linestrings] = + make_linestring_column({0, 2}, {0.0, 0.0, 1.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + CUSPATIAL_RUN_TEST(this->run_single, + geometry_column_view(linestrings->view(), ptype, geometry_type_id::LINESTRING), + geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON), + {1.0}); +}; + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, SingleToMultiOnePair) +{ + using T = TypeParam; + + auto [ptype, linestrings] = + make_linestring_column({0, 2}, {0.0, 0.0, 1.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + CUSPATIAL_RUN_TEST(this->run_single, + geometry_column_view(linestrings->view(), ptype, geometry_type_id::LINESTRING), + geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON), + {1.0}); +}; + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, MultiToSingleOnePair) +{ + using T = TypeParam; + + auto [ptype, linestrings] = + make_linestring_column({0, 1}, {0, 2}, {0.0, 0.0, 1.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + CUSPATIAL_RUN_TEST(this->run_single, + geometry_column_view(linestrings->view(), ptype, geometry_type_id::LINESTRING), + geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON), + {1.0}); +}; + +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, MultiToMultiOnePair) +{ + using T = TypeParam; + + auto [ptype, linestrings] = + make_linestring_column({0, 1}, {0, 2}, {0.0, 0.0, 1.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + CUSPATIAL_RUN_TEST(this->run_single, + geometry_column_view(linestrings->view(), ptype, geometry_type_id::LINESTRING), + geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON), + {1.0}); +}; + +TEST_F(PairwiseLinestringPolygonDistanceTestUntyped, SizeMismatch) +{ + auto [ptype, linestrings] = + make_linestring_column({0, 1, 2}, {0, 1, 2}, {0.0, 0.0, 1.0, 1.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto linestrings_view = + geometry_column_view(linestrings->view(), ptype, geometry_type_id::LINESTRING); + auto polygons_view = geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_linestring_polygon_distance(linestrings_view, polygons_view), + cuspatial::logic_error); +}; + +TEST_F(PairwiseLinestringPolygonDistanceTestUntyped, TypeMismatch) +{ + auto [ptype, linestrings] = + make_linestring_column({0, 1}, {0, 1}, {0.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto linestrings_view = + geometry_column_view(linestrings->view(), ptype, geometry_type_id::LINESTRING); + auto polygons_view = geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_linestring_polygon_distance(linestrings_view, polygons_view), + cuspatial::logic_error); +}; + +TEST_F(PairwiseLinestringPolygonDistanceTestUntyped, WrongGeometryType) +{ + auto [ptype, points] = make_point_column({0, 1}, {0.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto points_view = geometry_column_view(points->view(), ptype, geometry_type_id::POINT); + auto polygons_view = geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_linestring_polygon_distance(points_view, polygons_view), + cuspatial::logic_error); +}; diff --git a/cpp/tests/spatial/polygon_distance_test.cpp b/cpp/tests/spatial/polygon_distance_test.cpp new file mode 100644 index 000000000..a28d7a727 --- /dev/null +++ b/cpp/tests/spatial/polygon_distance_test.cpp @@ -0,0 +1,197 @@ +/* + * 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 +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include +#include + +using namespace cuspatial; +using namespace cuspatial::test; + +using namespace cudf; +using namespace cudf::test; + +template +struct PairwisePolygonDistanceTestBase : ::testing::Test { + void run_single(geometry_column_view lhs, + geometry_column_view rhs, + std::initializer_list expected) + { + auto got = pairwise_polygon_distance(lhs, rhs); + CUDF_TEST_EXPECT_COLUMNS_EQUIVALENT(*got, fixed_width_column_wrapper(expected)); + } + rmm::cuda_stream_view stream() { return cudf::get_default_stream(); } +}; + +template +struct PairwisePolygonDistanceTestEmpty : PairwisePolygonDistanceTestBase { + void SetUp() + { + [[maybe_unused]] collection_type_id _; + std::tie(_, empty_polygon_column) = make_polygon_column({0}, {0}, {}, this->stream()); + std::tie(_, empty_multipolygon_column) = + make_polygon_column({0}, {0}, {0}, {}, this->stream()); + } + + geometry_column_view empty_polygon() + { + return geometry_column_view( + empty_polygon_column->view(), collection_type_id::SINGLE, geometry_type_id::POLYGON); + } + + geometry_column_view empty_multipolygon() + { + return geometry_column_view( + empty_multipolygon_column->view(), collection_type_id::MULTI, geometry_type_id::POLYGON); + } + + std::unique_ptr empty_polygon_column; + std::unique_ptr empty_multipolygon_column; +}; + +using TestTypes = ::testing::Types; +TYPED_TEST_CASE(PairwisePolygonDistanceTestEmpty, TestTypes); +TYPED_TEST_CASE(PairwisePolygonDistanceTestOnePair, TestTypes); + +template +struct PairwisePolygonDistanceTestOnePair : PairwisePolygonDistanceTestBase { + void SetUp() + { + [[maybe_unused]] collection_type_id _; + std::tie(_, one_polygon_column) = + make_polygon_column({0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + std::tie(_, one_multipolygon_column) = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + } + + geometry_column_view one_polygon() + { + return geometry_column_view( + one_polygon_column->view(), collection_type_id::SINGLE, geometry_type_id::POLYGON); + } + + geometry_column_view one_multipolygon() + { + return geometry_column_view( + one_multipolygon_column->view(), collection_type_id::MULTI, geometry_type_id::POLYGON); + } + + std::unique_ptr one_polygon_column; + std::unique_ptr one_multipolygon_column; +}; + +struct PairwisePolygonDistanceTestUntyped : testing::Test { + rmm::cuda_stream_view stream() { return cudf::get_default_stream(); } +}; +TYPED_TEST(PairwisePolygonDistanceTestEmpty, SingleToSingleEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_polygon(), this->empty_polygon(), {}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestEmpty, SingleToMultiEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_polygon(), this->empty_multipolygon(), {}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestEmpty, MultiToSingleEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_multipolygon(), this->empty_polygon(), {}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestEmpty, MultiToMultiEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_multipolygon(), this->empty_multipolygon(), {}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestOnePair, SingleToSingleOnePair) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->one_polygon(), this->one_polygon(), {0}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestOnePair, SingleToMultiOnePair) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->one_polygon(), this->one_multipolygon(), {0}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestOnePair, MultiToSingleOnePair) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->one_multipolygon(), this->one_polygon(), {0}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestOnePair, MultiToMultiOnePair) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->one_multipolygon(), this->one_multipolygon(), {0}); +}; + +TEST_F(PairwisePolygonDistanceTestUntyped, SizeMismatch) +{ + auto [ptype, polygons1] = make_polygon_column( + {0, 1, 2}, {0, 4, 8}, {0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0}, this->stream()); + + auto [polytype, polygons2] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto polygons1_view = + geometry_column_view(polygons1->view(), ptype, geometry_type_id::LINESTRING); + auto polygons2_view = + geometry_column_view(polygons1->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_polygon_distance(polygons1_view, polygons2_view), cuspatial::logic_error); +}; + +TEST_F(PairwisePolygonDistanceTestUntyped, TypeMismatch) +{ + auto [ptype, polygons1] = make_polygon_column( + {0, 1, 2}, {0, 4, 8}, {0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0}, this->stream()); + + auto [polytype, polygons2] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto polygons1_view = + geometry_column_view(polygons1->view(), ptype, geometry_type_id::LINESTRING); + auto polygons2_view = + geometry_column_view(polygons2->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_polygon_distance(polygons1_view, polygons2_view), cuspatial::logic_error); +}; + +TEST_F(PairwisePolygonDistanceTestUntyped, WrongGeometryType) +{ + auto [ptype, points] = make_point_column({0, 1}, {0.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto points_view = geometry_column_view(points->view(), ptype, geometry_type_id::POINT); + auto polygons_view = geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_polygon_distance(points_view, polygons_view), cuspatial::logic_error); +}; diff --git a/python/cuspatial/cuspatial/__init__.py b/python/cuspatial/cuspatial/__init__.py index 2f3a27267..c72dc00bb 100644 --- a/python/cuspatial/cuspatial/__init__.py +++ b/python/cuspatial/cuspatial/__init__.py @@ -7,10 +7,12 @@ join_quadtree_and_bounding_boxes, linestring_bounding_boxes, pairwise_linestring_distance, + pairwise_linestring_polygon_distance, pairwise_point_distance, pairwise_point_linestring_distance, pairwise_point_linestring_nearest_points, pairwise_point_polygon_distance, + pairwise_polygon_distance, point_in_polygon, points_in_spatial_window, polygon_bounding_boxes, diff --git a/python/cuspatial/cuspatial/_lib/cpp/distance/linestring_polygon_distance.pxd b/python/cuspatial/cuspatial/_lib/cpp/distance/linestring_polygon_distance.pxd new file mode 100644 index 000000000..82a4f1834 --- /dev/null +++ b/python/cuspatial/cuspatial/_lib/cpp/distance/linestring_polygon_distance.pxd @@ -0,0 +1,17 @@ +# Copyright (c) 2023, NVIDIA CORPORATION. + +from libcpp.memory cimport unique_ptr + +from cudf._lib.cpp.column.column cimport column + +from cuspatial._lib.cpp.column.geometry_column_view cimport ( + geometry_column_view, +) + + +cdef extern from "cuspatial/distance/linestring_polygon_distance.hpp" \ + namespace "cuspatial" nogil: + cdef unique_ptr[column] pairwise_linestring_polygon_distance( + const geometry_column_view & multilinestrings, + const geometry_column_view & multipolygons + ) except + diff --git a/python/cuspatial/cuspatial/_lib/cpp/distance/polygon_distance.pxd b/python/cuspatial/cuspatial/_lib/cpp/distance/polygon_distance.pxd new file mode 100644 index 000000000..9d2bfab39 --- /dev/null +++ b/python/cuspatial/cuspatial/_lib/cpp/distance/polygon_distance.pxd @@ -0,0 +1,17 @@ +# Copyright (c) 2023, NVIDIA CORPORATION. + +from libcpp.memory cimport unique_ptr + +from cudf._lib.cpp.column.column cimport column + +from cuspatial._lib.cpp.column.geometry_column_view cimport ( + geometry_column_view, +) + + +cdef extern from "cuspatial/distance/polygon_distance.hpp" \ + namespace "cuspatial" nogil: + cdef unique_ptr[column] pairwise_polygon_distance( + const geometry_column_view & lhs, + const geometry_column_view & rhs + ) except + diff --git a/python/cuspatial/cuspatial/_lib/distance.pyx b/python/cuspatial/cuspatial/_lib/distance.pyx index f3d68717e..7edf19754 100644 --- a/python/cuspatial/cuspatial/_lib/distance.pyx +++ b/python/cuspatial/cuspatial/_lib/distance.pyx @@ -13,6 +13,9 @@ from cuspatial._lib.cpp.column.geometry_column_view cimport ( from cuspatial._lib.cpp.distance.linestring_distance cimport ( pairwise_linestring_distance as c_pairwise_linestring_distance, ) +from cuspatial._lib.cpp.distance.linestring_polygon_distance cimport ( + pairwise_linestring_polygon_distance as c_pairwise_line_poly_dist, +) from cuspatial._lib.cpp.distance.point_distance cimport ( pairwise_point_distance as c_pairwise_point_distance, ) @@ -22,6 +25,9 @@ from cuspatial._lib.cpp.distance.point_linestring_distance cimport ( from cuspatial._lib.cpp.distance.point_polygon_distance cimport ( pairwise_point_polygon_distance as c_pairwise_point_polygon_distance, ) +from cuspatial._lib.cpp.distance.polygon_distance cimport ( + pairwise_polygon_distance as c_pairwise_polygon_distance, +) from cuspatial._lib.cpp.optional cimport optional from cuspatial._lib.cpp.types cimport collection_type_id, geometry_type_id from cuspatial._lib.types cimport collection_type_py_to_c @@ -142,3 +148,53 @@ def pairwise_point_polygon_distance( )) return Column.from_unique_ptr(move(c_result)) + + +def pairwise_linestring_polygon_distance( + Column multilinestrings, + Column multipolygons +): + + cdef shared_ptr[geometry_column_view] c_multilinestrings = \ + make_shared[geometry_column_view]( + multilinestrings.view(), + collection_type_id.MULTI, + geometry_type_id.LINESTRING) + + cdef shared_ptr[geometry_column_view] c_multipolygons = \ + make_shared[geometry_column_view]( + multipolygons.view(), + collection_type_id.MULTI, + geometry_type_id.POLYGON) + + cdef unique_ptr[column] c_result + + with nogil: + c_result = move(c_pairwise_line_poly_dist( + c_multilinestrings.get()[0], c_multipolygons.get()[0] + )) + + return Column.from_unique_ptr(move(c_result)) + + +def pairwise_polygon_distance(Column lhs, Column rhs): + cdef shared_ptr[geometry_column_view] c_lhs = \ + make_shared[geometry_column_view]( + lhs.view(), + collection_type_id.MULTI, + geometry_type_id.POLYGON) + + cdef shared_ptr[geometry_column_view] c_rhs = \ + make_shared[geometry_column_view]( + rhs.view(), + collection_type_id.MULTI, + geometry_type_id.POLYGON) + + cdef unique_ptr[column] c_result + + with nogil: + c_result = move(c_pairwise_polygon_distance( + c_lhs.get()[0], c_rhs.get()[0] + )) + + return Column.from_unique_ptr(move(c_result)) diff --git a/python/cuspatial/cuspatial/core/binops/distance.py b/python/cuspatial/cuspatial/core/binops/distance.py new file mode 100644 index 000000000..fcd868190 --- /dev/null +++ b/python/cuspatial/cuspatial/core/binops/distance.py @@ -0,0 +1,47 @@ +# Copyright (c) 2023, NVIDIA CORPORATION. + +from itertools import product + +from cuspatial._lib.distance import ( + pairwise_point_distance , + pairwise_point_linestring_distance , + pairwise_point_polygon_distance , + pairwise_linestring_distance , + pairwise_linestring_polygon_distance , + pairwise_polygon_distance , +) + +from cuspatial.core._column.geometa import Feature_Enum as F +from cuspatial.core._column.geocolumn import GeoColumn +from cuspatial.core.dispatch.geocolumn_binop_dispatch import GeoColumnBinopDispatch + + +class DistanceDispatch(GeoColumnBinopDispatch): + def __init__(self, lhs: GeoColumn, rhs: GeoColumn): + super().__init__(lhs, rhs) + self.dispatch_dict = { + (F.POINT, F.POINT): (pairwise_point_distance, False), + (F.POINT, F.MULTIPOINT): (pairwise_point_distance, False), + (F.POINT, F.LINESTRING): (pairwise_point_linestring_distance, False), + (F.POINT, F.POLYGON): (pairwise_point_polygon_distance, False), + + (F.MULTIPOINT, F.POINT): (pairwise_point_distance, False), + (F.MULTIPOINT, F.MULTIPOINT): (pairwise_point_distance, False), + (F.MULTIPOINT, F.LINESTRING): (pairwise_point_linestring_distance, False), + (F.MULTIPOINT, F.POLYGON): (pairwise_point_polygon_distance, False), + + (F.LINESTRING, F.POINT): (pairwise_point_linestring_distance, True), + (F.LINESTRING, F.MULTIPOINT): (pairwise_point_linestring_distance, True), + (F.LINESTRING, F.LINESTRING): (pairwise_linestring_distance, False), + (F.LINESTRING, F.POLYGON): (pairwise_linestring_polygon_distance, False), + + (F.POLYGON, F.POINT): (pairwise_point_polygon_distance, True), + (F.POLYGON, F.MULTIPOINT): (pairwise_point_polygon_distance, True), + (F.POLYGON, F.LINESTRING): (pairwise_linestring_polygon_distance, True), + (F.POLYGON, F.POLYGON): (pairwise_polygon_distance, False), + } + + none = [F.NONE] + other = [F.POINT, F.MULTIPOINT, F.LINESTRING, F.POLYGON] + for ltype, rtype in [*product(none, other), *product(other, none)]: + self.dispatch_dict[(ltype, rtype)] = ("Impossible", None) diff --git a/python/cuspatial/cuspatial/core/dispatch/geocolumn_binop_dispatch.py b/python/cuspatial/cuspatial/core/dispatch/geocolumn_binop_dispatch.py new file mode 100644 index 000000000..aa18e0cf0 --- /dev/null +++ b/python/cuspatial/cuspatial/core/dispatch/geocolumn_binop_dispatch.py @@ -0,0 +1,118 @@ +import cupy as cp + +import cudf + +from cuspatial.core._column.geocolumn import GeoColumn +from cuspatial.core._column.geometa import Feature_Enum + + +class GeoColumnBinopDispatch: + """Dispatch binary operations of two geocolumns to subkernels and + reassemble data. + """ + + def __init__(self, lhs: GeoColumn, rhs: GeoColumn): + if len(lhs) != len(rhs): + raise ValueError("Input is of different length.") + + self.lhs = lhs + self.rhs = rhs + self.dispatch_dict = None + + def __call__(self): + if self.dispatch_dict is None: + raise NotImplementedError( + "Binary op dispatch base class cannot be used directly, use a " + "derived class instead." + ) + self._sort_geocolumn_by_type_pairs() + self._compute_gather_offsets() + self._gather_and_dispatch_computation() + self._assemble_results() + return self._reordered_result + + def _sort_geocolumn_by_type_pairs(self): + """Given two geocolumn lhs and rhs, sort the offset buffer using + the pair of types buffer column as keys. + + This result in a reordered view into child column where each subtype + is contiguous in the offset buffer. + """ + + df = cudf.DataFrame( + { + "lhs_types": self.lhs._meta.input_types, + "rhs_types": self.rhs._meta.input_types, + "lhs_offsets": self.lhs._meta.union_offsets, + "rhs_offsets": self.rhs._meta.union_offsets, + "order": cp.arange(len(self.lhs)), + } + ) + + self.df_sorted = df.sort_values(by=["lhs_types", "rhs_types"]) + + def _compute_gather_offsets(self): + """From the sorted types buffer and offset buffer, compute the + boundaries of each contiguous section. + """ + self.binary_type_to_offsets = {} + self.orders = [] + for lhs_type, rhs_type in self.dispatch_dict: + mask = ( + (self.df_sorted["lhs_types"] + == lhs_type.value) & (self.df_sorted["rhs_types"] + == rhs_type.value) + ) + masked_df = self.df_sorted[["lhs_offsets", "rhs_offsets", "order"]][ + mask + ] + self.binary_type_to_offsets[(lhs_type, rhs_type)] = ( + masked_df["lhs_offsets"]._column, + masked_df["rhs_offsets"]._column, + ) + self.orders.append(masked_df["order"]) + + def _gather_and_dispatch_computation(self): + """For each type combination, gather from the child columns and + dispatch to child kernel. + """ + self.results = [] + for lhs_type, rhs_type in self.dispatch_dict: + op, reflect = self.dispatch_dict[(lhs_type, rhs_type)] + lhs_offsets, rhs_offsets = self.binary_type_to_offsets[ + (lhs_type, rhs_type) + ] + if op == "Impossible": + self.result[(lhs_type, rhs_type)] = cp.full( + (len(lhs_offsets)), float("nan"), "f8" + ) + else: + lhs = self._gather_by_type(self.lhs, lhs_type, lhs_offsets) + rhs = self._gather_by_type(self.rhs, rhs_type, rhs_offsets) + self.results.append( + op(lhs, rhs) if not reflect else op(rhs, lhs) + ) + + def _assemble_results(self): + """After computed all results from each subkernel, concatenate the + result and make sure that the order matches that of original input. + """ + if len(self.results) == 0: + self.reordered_results = cudf.Series() + + results_concat = cudf.concat(self.results) + orders_concat = cudf.concat(self.orders) + self.reordered_results = cudf.Series(cp.zeros(results_concat)) + self.reordered_results[orders_concat] = results_concat + + def _gather_by_type(self, column, typ, offsets): + if typ == Feature_Enum.POINT: + return column.points._column.take(offsets) + elif typ == Feature_Enum.MULTIPOINT: + return column.mpoints._column.take(offsets) + elif typ == Feature_Enum.LINESTRING: + return column.lines._column.take(offsets) + elif typ == Feature_Enum.POLYGON: + return column.polygons._column.take(offsets) + else: + raise ValueError(f"Unrecognized type enum: {typ}") diff --git a/python/cuspatial/cuspatial/core/geoseries.py b/python/cuspatial/cuspatial/core/geoseries.py index 9ff85a7dc..50d37cc6b 100644 --- a/python/cuspatial/cuspatial/core/geoseries.py +++ b/python/cuspatial/cuspatial/core/geoseries.py @@ -42,6 +42,7 @@ contains_only_points, contains_only_polygons, ) +from cuspatial.core.binops.distance import DistanceDispatch T = TypeVar("T", bound="GeoSeries") @@ -1240,3 +1241,29 @@ def disjoint(self, other, align=True): align=align ) return predicate(self, other) + + + def distance(self, other, align=True): + """Returns a Series containing the distance to aligned other. + + Parameters + ---------- + Parameters + other: cuspatial.GeoSeries + The Geoseries (elementwise) or geometric object to find the + distance to. + + align: bool (default True) + If True, automatically aligns GeoSeries based on their indices. + If False, the order of elements is preserved. + + Returns + ------- + result : cudf.Series + A series of floating point values of the distance between the + corresponding pair + """ + if align: + lhs, rhs = self.align(other) + + return DistanceDispatch(lhs._column, rhs._column)() diff --git a/python/cuspatial/cuspatial/core/spatial/__init__.py b/python/cuspatial/cuspatial/core/spatial/__init__.py index ee07838f9..97486a877 100644 --- a/python/cuspatial/cuspatial/core/spatial/__init__.py +++ b/python/cuspatial/cuspatial/core/spatial/__init__.py @@ -5,9 +5,11 @@ directed_hausdorff_distance, haversine_distance, pairwise_linestring_distance, + pairwise_linestring_polygon_distance, pairwise_point_distance, pairwise_point_linestring_distance, pairwise_point_polygon_distance, + pairwise_polygon_distance, ) from .filtering import points_in_spatial_window from .indexing import quadtree_on_points @@ -27,9 +29,11 @@ "sinusoidal_projection", "pairwise_point_distance", "pairwise_linestring_distance", + "pairwise_linestring_polygon_distance", "pairwise_point_polygon_distance", "pairwise_point_linestring_distance", "pairwise_point_linestring_nearest_points", + "pairwise_polygon_distance", "polygon_bounding_boxes", "linestring_bounding_boxes", "point_in_polygon", diff --git a/python/cuspatial/cuspatial/core/spatial/distance.py b/python/cuspatial/cuspatial/core/spatial/distance.py index f48f6a432..2e7ff976f 100644 --- a/python/cuspatial/cuspatial/core/spatial/distance.py +++ b/python/cuspatial/cuspatial/core/spatial/distance.py @@ -8,9 +8,11 @@ from cuspatial._lib.distance import ( pairwise_linestring_distance as cpp_pairwise_linestring_distance, + pairwise_linestring_polygon_distance as c_pairwise_line_poly_dist, pairwise_point_distance as cpp_pairwise_point_distance, pairwise_point_linestring_distance as c_pairwise_point_linestring_distance, pairwise_point_polygon_distance as c_pairwise_point_polygon_distance, + pairwise_polygon_distance as c_pairwise_polygon_distance, ) from cuspatial._lib.hausdorff import ( directed_hausdorff_distance as cpp_directed_hausdorff_distance, @@ -443,7 +445,7 @@ def pairwise_point_polygon_distance(points: GeoSeries, polygons: GeoSeries): raise ValueError("`points` array must contain only points") if not contains_only_polygons(polygons): - raise ValueError("`linestrings` array must contain only linestrings") + raise ValueError("`polygons` array must contain only polygons") if len(points.points.xy) > 0 and len(points.multipoints.xy) > 0: raise NotImplementedError( @@ -480,6 +482,149 @@ def pairwise_point_polygon_distance(points: GeoSeries, polygons: GeoSeries): ) +def pairwise_linestring_polygon_distance( + linestrings: GeoSeries, polygons: GeoSeries +): + """Compute distance between pairs of (multi)linestrings and (multi)polygons + + The distance between a (multi)linestrings and a (multi)polygon + is defined as the shortest distance between every segment in the + multilinestring and every edge of the (multi)polygon. If the + multilinestring and multipolygon intersects, the distance is 0. + + This algorithm computes distance pairwise. The ith row in the result is + the distance between the ith (multi)linestring in `linestrings` and the ith + (multi)polygon in `polygons`. + + Parameters + ---------- + linestrings : GeoSeries + The (multi)linestrings to compute the distance from. + polygons : GeoSeries + The (multi)polygons to compute the distance from. + + Returns + ------- + distance : cudf.Series + + Notes + ----- + The input `GeoSeries` must contain a single type geometry. + For example, `linestrings` series cannot contain both linestrings and + polygons. + + Examples + -------- + Compute distance between a linestring and a polygon: + >>> from shapely.geometry import LineString, Polygon + >>> lines = cuspatial.GeoSeries([ + ... LineString([(0, 0), (1, 1)])]) + >>> polys = cuspatial.GeoSeries([ + ... Polygon([(-1, -1), (-1, 0), (-2, 0), (-1, -1)]) + ... ]) + >>> cuspatial.pairwise_linestring_polygon_distance(lines, polys) + 0 1.0 + dtype: float64 + + Compute distance between a multipoint and a multipolygon + >>> from shapely.geometry import MultiLineString, MultiPolygon + >>> lines = cuspatial.GeoSeries([ + ... MultiLineString([ + ... LineString([(0, 0), (1, 1)]), + ... LineString([(1, 1), (2, 2)])]) + ... ]) + >>> polys = cuspatial.GeoSeries([ + ... MultiPolygon([ + ... Polygon([(-1, -1), (-1, 0), (-2, 0), (-1, -1)]), + ... Polygon([(-2, 0), (-3, 0), (-3, -1), (-2, 0)])]) + ... ]) + >>> cuspatial.pairwise_linestring_polygon_distance(lines, polys) + 0 1.0 + dtype: float64 + """ + + if len(linestrings) != len(polygons): + raise ValueError("Unmatched input geoseries length.") + + if len(linestrings) == 0: + return cudf.Series(dtype=linestrings.lines.xy.dtype) + + if not contains_only_linestrings(linestrings): + raise ValueError("`linestrings` array must contain only linestrings") + + if not contains_only_polygons(polygons): + raise ValueError("`polygon` array must contain only polygons") + + # Handle slicing in geoseries + linestrings_column = linestrings._column.lines._column.take( + linestrings._column._meta.union_offsets._column + ) + polygon_column = polygons._column.polygons._column.take( + polygons._column._meta.union_offsets._column + ) + + return Series._from_data( + {None: c_pairwise_line_poly_dist(linestrings_column, polygon_column)} + ) + + +def pairwise_polygon_distance(polygons1: GeoSeries, polygons2: GeoSeries): + """Compute distance between pairs of (multi)polygons and (multi)polygons + The distance between a (multi)polygon and a (multi)polygon + is defined as the shortest distance between every edge of the + (multi)polygon pair. If the multipolygon and multipolygon intersects, + the distance is 0. + + This algorithm computes distance pairwise. The ith row in the result is + the distance between the ith (multi)polygon in `polygons1` and the ith + (multi)polygon in `polygons2`. + + Parameters + ---------- + polygons1 : GeoSeries + The (multi)polygons to compute the distance from. + polygons2 : GeoSeries + The (multi)polygons to compute the distance from. + Returns + ------- + distance : cudf.Series + + Notes + ----- + The input `GeoSeries` must contain a single type geometry. + For example, `polygons1` series cannot contain both points and + polygons. + + Examples + -------- + Compute distance between a point and a polygon: + """ + + if len(polygons1) != len(polygons2): + raise ValueError("Unmatched input geoseries length.") + + if len(polygons1) == 0: + return cudf.Series(dtype=polygons1.lines.xy.dtype) + + if not contains_only_polygons(polygons1): + raise ValueError("`polygons1` array must contain only polygons") + + if not contains_only_polygons(polygons2): + raise ValueError("`polygons2` array must contain only polygons") + + # Handle slicing and aligns in geoseries + polygon1_column = polygons1._column.polygons._column.take( + polygons1._column._meta.union_offsets._column + ) + polygon2_column = polygons2._column.polygons._column.take( + polygons2._column._meta.union_offsets._column + ) + + return Series._from_data( + {None: c_pairwise_polygon_distance(polygon1_column, polygon2_column)} + ) + + def _flatten_point_series( points: GeoSeries, ) -> Tuple[ diff --git a/python/cuspatial/cuspatial/tests/spatial/distance/test_distance.py b/python/cuspatial/cuspatial/tests/spatial/distance/test_distance.py new file mode 100644 index 000000000..be5053320 --- /dev/null +++ b/python/cuspatial/cuspatial/tests/spatial/distance/test_distance.py @@ -0,0 +1,57 @@ +import pytest + +from shapely.geometry import * + +from cudf.testing import assert_series_equal + +import cuspatial + +@pytest.mark.parametrize( + "lobj", [ + Point(0, 0), + MultiPoint([ + Point(1, 1), + Point(0, 0) + ]), + LineString([(0, 0), (1, 1)]), + MultiLineString([ + LineString([(0, 0), (1, 1)]), + LineString([(2, 2), (3, 3)]), + ]), + Polygon([(0, 0), (1, 0), (0, 1), (0, 0)]), + MultiPolygon([ + Polygon([(0, 0), (1, 0), (0, 1), (0, 0)]), + Polygon([(2, 2), (2, 3), (3, 2), (2, 2)]) + ]) + ] +) +@pytest.mark.parametrize( + "robj", [ + Point(10, 10), + MultiPoint([ + Point(11, 11), + Point(10, 10) + ]), + LineString([(10, 10), (11, 11)]), + MultiLineString([ + LineString([(10, 10), (11, 11)]), + LineString([(12, 12), (13, 13)]), + ]), + Polygon([(10, 0), (11, 10), (10, 11), (10, 10)]), + MultiPolygon([ + Polygon([(10, 10), (11, 10), (10, 11), (10, 10)]), + Polygon([(12, 12), (12, 13), (13, 12), (12, 12)]) + ]) + ] +) +def test_one_pair_all_combos_no_nulls_test(lobj, robj): + lhs = cuspatial.GeoSeries([lobj]) + rhs = cuspatial.GeoSeries([robj]) + + hlhs = lhs.to_geopandas() + hrhs = rhs.to_geopandas() + + got = lhs.distance(rhs) + expect = hlhs.distance(hrhs) + + assert_series_equal(cudf.Series(expect), got) diff --git a/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_linestring_polygon_distance.py b/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_linestring_polygon_distance.py new file mode 100644 index 000000000..f092b6023 --- /dev/null +++ b/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_linestring_polygon_distance.py @@ -0,0 +1,117 @@ +import geopandas as gpd +import pytest +from shapely.geometry import LineString, MultiLineString, MultiPolygon, Polygon + +import cudf +from cudf.testing import assert_series_equal + +import cuspatial + + +def test_linestring_polygon_empty(): + lhs = cuspatial.GeoSeries.from_linestrings_xy([], [0], [0]) + rhs = cuspatial.GeoSeries.from_polygons_xy([], [0], [0], [0]) + + got = cuspatial.pairwise_linestring_polygon_distance(lhs, rhs) + + expect = cudf.Series([], dtype="f8") + + assert_series_equal(got, expect) + + +@pytest.mark.parametrize( + "linestrings", + [ + [LineString([(0, 0), (1, 1)])], + [MultiLineString([[(1, 1), (2, 2)], [(10, 10), (11, 11)]])], + ], +) +@pytest.mark.parametrize( + "polygons", + [ + [Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)])], + [ + MultiPolygon( + [ + Polygon([(-2, 0), (-1, 0), (-1, -1), (-2, 0)]), + Polygon([(1, 0), (2, 0), (1, -1), (1, 0)]), + ] + ) + ], + ], +) +def test_one_pair(linestrings, polygons): + lhs = gpd.GeoSeries(linestrings) + rhs = gpd.GeoSeries(polygons) + + dlhs = cuspatial.GeoSeries(linestrings) + drhs = cuspatial.GeoSeries(polygons) + + expect = lhs.distance(rhs) + got = cuspatial.pairwise_linestring_polygon_distance(dlhs, drhs) + + assert_series_equal(got, cudf.Series(expect)) + + +@pytest.mark.parametrize( + "linestrings", + [ + [LineString([(0, 0), (1, 1)]), LineString([(10, 10), (11, 11)])], + [ + MultiLineString([[(1, 1), (2, 2)], [(3, 3), (4, 4)]]), + MultiLineString([[(10, 10), (11, 11)], [(12, 12), (13, 13)]]), + ], + ], +) +@pytest.mark.parametrize( + "polygons", + [ + [ + Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)]), + Polygon([(-4, -4), (-4, -5), (-5, -5), (-5, -4), (-5, -5)]), + ], + [ + MultiPolygon( + [ + Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)]), + Polygon([(0, 1), (1, 0), (0, -1), (-1, 0), (0, 1)]), + ] + ), + MultiPolygon( + [ + Polygon( + [(-4, -4), (-4, -5), (-5, -5), (-5, -4), (-5, -5)] + ), + Polygon([(-2, 0), (-2, -2), (0, -2), (0, 0), (-2, 0)]), + ] + ), + ], + ], +) +def test_two_pair(linestrings, polygons): + lhs = gpd.GeoSeries(linestrings) + rhs = gpd.GeoSeries(polygons) + + dlhs = cuspatial.GeoSeries(linestrings) + drhs = cuspatial.GeoSeries(polygons) + + expect = lhs.distance(rhs) + got = cuspatial.pairwise_linestring_polygon_distance(dlhs, drhs) + + assert_series_equal(got, cudf.Series(expect)) + + +def test_linestring_polygon_large(linestring_generator, polygon_generator): + N = 100 + linestrings = gpd.GeoSeries(linestring_generator(N, 5)) + polygons = gpd.GeoSeries(polygon_generator(N, 10.0, 3.0)) + + dlinestrings = cuspatial.from_geopandas(linestrings) + dpolygons = cuspatial.from_geopandas(polygons) + + expect = linestrings.distance(polygons) + got = cuspatial.pairwise_linestring_polygon_distance( + dlinestrings, dpolygons + ) + + assert_series_equal(got, cudf.Series(expect)) diff --git a/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_point_polygon_distance.py b/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_point_polygon_distance.py index 199c41208..49fabec39 100644 --- a/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_point_polygon_distance.py +++ b/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_point_polygon_distance.py @@ -124,13 +124,6 @@ def test_point_polygon_geocities(naturalearth_cities, naturalearth_lowres): gpu_cities = cuspatial.from_geopandas(naturalearth_cities.geometry) gpu_countries = cuspatial.from_geopandas(naturalearth_lowres.geometry) - print( - len(naturalearth_lowres), - len(naturalearth_lowres[: len(naturalearth_cities)]), - len(gpu_countries), - len(gpu_countries[: len(naturalearth_cities)]), - ) - expect = naturalearth_cities.geometry[:N].distance( naturalearth_lowres.geometry[:N] ) diff --git a/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_polygon_distance.py b/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_polygon_distance.py new file mode 100644 index 000000000..12095bebd --- /dev/null +++ b/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_polygon_distance.py @@ -0,0 +1,149 @@ +import geopandas as gpd +import pytest +from shapely.geometry import MultiPolygon, Polygon + +import cudf +from cudf.testing import assert_series_equal + +import cuspatial + + +def test_polygon_empty(): + lhs = cuspatial.GeoSeries.from_polygons_xy([], [0], [0], [0]) + rhs = cuspatial.GeoSeries.from_polygons_xy([], [0], [0], [0]) + + got = cuspatial.pairwise_polygon_distance(lhs, rhs) + + expect = cudf.Series([], dtype="f8") + + assert_series_equal(got, expect) + + +@pytest.mark.parametrize( + "polygons1", + [ + [Polygon([(10, 11), (11, 10), (11, 11), (10, 11)])], + [ + MultiPolygon( + [ + Polygon([(12, 10), (11, 10), (11, 11), (12, 10)]), + Polygon([(11, 10), (12, 10), (11, 11), (11, 10)]), + ] + ) + ], + ], +) +@pytest.mark.parametrize( + "polygons2", + [ + [Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)])], + [ + MultiPolygon( + [ + Polygon([(-2, 0), (-1, 0), (-1, -1), (-2, 0)]), + Polygon([(1, 0), (2, 0), (1, -1), (1, 0)]), + ] + ) + ], + ], +) +def test_one_pair(polygons1, polygons2): + lhs = gpd.GeoSeries(polygons1) + rhs = gpd.GeoSeries(polygons2) + + dlhs = cuspatial.GeoSeries(polygons1) + drhs = cuspatial.GeoSeries(polygons2) + + expect = lhs.distance(rhs) + got = cuspatial.pairwise_polygon_distance(dlhs, drhs) + + assert_series_equal(got, cudf.Series(expect)) + + +@pytest.mark.parametrize( + "polygons1", + [ + [ + Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)]), + Polygon([(-4, -4), (-4, -5), (-5, -5), (-5, -4), (-5, -5)]), + ], + [ + MultiPolygon( + [ + Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)]), + Polygon([(0, 1), (1, 0), (0, -1), (-1, 0), (0, 1)]), + ] + ), + MultiPolygon( + [ + Polygon( + [(-4, -4), (-4, -5), (-5, -5), (-5, -4), (-5, -5)] + ), + Polygon([(-2, 0), (-2, -2), (0, -2), (0, 0), (-2, 0)]), + ] + ), + ], + ], +) +@pytest.mark.parametrize( + "polygons2", + [ + [ + Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)]), + Polygon([(-4, -4), (-4, -5), (-5, -5), (-5, -4), (-5, -5)]), + ], + [ + MultiPolygon( + [ + Polygon([(0, 1), (1, 0), (-1, 0), (0, 1)]), + Polygon([(0, 1), (1, 0), (0, -1), (-1, 0), (0, 1)]), + ] + ), + MultiPolygon( + [ + Polygon( + [(-4, -4), (-4, -5), (-5, -5), (-5, -4), (-5, -5)] + ), + Polygon([(-2, 0), (-2, -2), (0, -2), (0, 0), (-2, 0)]), + ] + ), + ], + ], +) +def test_two_pair(polygons1, polygons2): + lhs = gpd.GeoSeries(polygons1) + rhs = gpd.GeoSeries(polygons2) + + dlhs = cuspatial.GeoSeries(polygons1) + drhs = cuspatial.GeoSeries(polygons2) + + expect = lhs.distance(rhs) + got = cuspatial.pairwise_polygon_distance(dlhs, drhs) + + assert_series_equal(got, cudf.Series(expect)) + + +def test_linestring_polygon_large(polygon_generator): + N = 100 + polygons1 = gpd.GeoSeries(polygon_generator(N, 20.0, 5.0)) + polygons2 = gpd.GeoSeries(polygon_generator(N, 10.0, 3.0)) + + dpolygons1 = cuspatial.from_geopandas(polygons1) + dpolygons2 = cuspatial.from_geopandas(polygons2) + + expect = polygons1.distance(polygons2) + got = cuspatial.pairwise_polygon_distance(dpolygons1, dpolygons2) + + assert_series_equal(got, cudf.Series(expect)) + + +def test_point_polygon_geoboundaries(naturalearth_lowres): + N = 50 + + lhs = naturalearth_lowres.geometry[:N].reset_index(drop=True) + rhs = naturalearth_lowres.geometry[N : 2 * N].reset_index(drop=True) + expect = lhs.distance(rhs) + got = cuspatial.pairwise_polygon_distance( + cuspatial.GeoSeries(lhs), cuspatial.GeoSeries(rhs) + ) + assert_series_equal(cudf.Series(expect), got)