diff --git a/docs/source/user_guide/cuspatial_api_examples.ipynb b/docs/source/user_guide/cuspatial_api_examples.ipynb index 761b11831..bd78b237f 100644 --- a/docs/source/user_guide/cuspatial_api_examples.ipynb +++ b/docs/source/user_guide/cuspatial_api_examples.ipynb @@ -889,6 +889,161 @@ "print(gpu_polygons.head())" ] }, + { + "cell_type": "markdown", + "id": "008d320d-ca47-459f-9fff-8769494c8a61", + "metadata": {}, + "source": [ + "### cuspatial.pairwise_point_polygon_distance\n", + "\n", + "Using WGS 84 Pseudo-Mercator, distances are in meters." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "258c9a8c-7fe3-4047-80b7-00878d9fb2f1", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
pop_estcontinentnameiso_a3gdp_md_estgeometrydistance_fromdistance
0889953.0OceaniaFijiFJI5496MULTIPOLYGON (((20037508.343 -1812498.413, 200...Vatican City1.969350e+07
158005463.0AfricaTanzaniaTZA63177POLYGON ((3774143.866 -105758.362, 3792946.708...San Marino5.929777e+06
2603253.0AfricaW. SaharaESH907POLYGON ((-964649.018 3205725.605, -964597.245...Vaduz3.421172e+06
337589262.0North AmericaCanadaCAN1736425MULTIPOLYGON (((-13674486.249 6274861.394, -13...Lobamba1.296059e+07
4328239523.0North AmericaUnited States of AmericaUSA21433226MULTIPOLYGON (((-13674486.249 6274861.394, -13...Luxembourg8.174897e+06
\n", + "
" + ], + "text/plain": [ + " pop_est continent name iso_a3 gdp_md_est \\\n", + "0 889953.0 Oceania Fiji FJI 5496 \n", + "1 58005463.0 Africa Tanzania TZA 63177 \n", + "2 603253.0 Africa W. Sahara ESH 907 \n", + "3 37589262.0 North America Canada CAN 1736425 \n", + "4 328239523.0 North America United States of America USA 21433226 \n", + "\n", + " geometry distance_from \\\n", + "0 MULTIPOLYGON (((20037508.343 -1812498.413, 200... Vatican City \n", + "1 POLYGON ((3774143.866 -105758.362, 3792946.708... San Marino \n", + "2 POLYGON ((-964649.018 3205725.605, -964597.245... Vaduz \n", + "3 MULTIPOLYGON (((-13674486.249 6274861.394, -13... Lobamba \n", + "4 MULTIPOLYGON (((-13674486.249 6274861.394, -13... Luxembourg \n", + "\n", + " distance \n", + "0 1.969350e+07 \n", + "1 5.929777e+06 \n", + "2 3.421172e+06 \n", + "3 1.296059e+07 \n", + "4 8.174897e+06 \n", + "(GPU)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cities = geopandas.read_file(geopandas.datasets.get_path(\"naturalearth_cities\")).to_crs(3857)\n", + "countries = geopandas.read_file(geopandas.datasets.get_path(\"naturalearth_lowres\")).to_crs(3857)\n", + "\n", + "gpu_cities = cuspatial.from_geopandas(cities)\n", + "gpu_countries = cuspatial.from_geopandas(countries)\n", + "\n", + "dist = cuspatial.pairwise_point_polygon_distance(\n", + " gpu_cities.geometry[:len(gpu_countries)], gpu_countries.geometry\n", + ")\n", + "\n", + "gpu_countries[\"distance_from\"] = cities.name\n", + "gpu_countries[\"distance\"] = dist\n", + "\n", + "gpu_countries.head()" + ] + }, { "attachments": { "351aea0c-f37e-4ab9-bad2-c67bce69b5c3.png": { @@ -910,7 +1065,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 18, "id": "d1ade9da-c9e2-45c4-9685-dffeda3fd358", "metadata": { "tags": [] @@ -975,7 +1130,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 19, "id": "cc72a44d-a9bf-4432-9898-de899ac45869", "metadata": { "tags": [] @@ -993,7 +1148,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 20, "id": "1125fd17-afe1-4b9c-b48d-8842dd3700b3", "metadata": { "tags": [] @@ -1002,7 +1157,7 @@ { "data": { "text/plain": [ - "\n", + "\n", "[\n", " 0,\n", " 144\n", @@ -1010,7 +1165,7 @@ "dtype: int32" ] }, - "execution_count": 19, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -1023,7 +1178,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 21, "id": "b281e3bb-42d4-4d60-9cb2-b7dcc20b4776", "metadata": { "tags": [] @@ -1046,7 +1201,7 @@ "Length: 144, dtype: geometry" ] }, - "execution_count": 20, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } @@ -1058,7 +1213,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 22, "id": "e19873b9-2614-4242-ad67-caa47f807d04", "metadata": { "tags": [] @@ -1117,7 +1272,7 @@ "0 [9, 10, 10, 11, 11, 28, 12, 12, 13, 13, 14, 15... " ] }, - "execution_count": 21, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } @@ -1150,7 +1305,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 23, "id": "bf7b2256", "metadata": { "tags": [] @@ -1167,7 +1322,7 @@ "dtype: int64" ] }, - "execution_count": 22, + "execution_count": 23, "metadata": {}, "output_type": "execute_result" } @@ -1252,7 +1407,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 24, "id": "e3a0a9a3-0bdd-4f05-bcb5-7db4b99a44a3", "metadata": { "tags": [] @@ -1316,7 +1471,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 25, "id": "023bd25a-35be-435d-ab0b-ecbd7a47e147", "metadata": { "tags": [] @@ -1375,7 +1530,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 26, "id": "784aff8e-c9ed-4a81-aa87-bf301b3b90af", "metadata": { "tags": [] @@ -1390,7 +1545,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 27, "id": "fea24c78-cf5c-45c6-b860-338238e61323", "metadata": { "tags": [] diff --git a/python/cuspatial/cuspatial/__init__.py b/python/cuspatial/cuspatial/__init__.py index bb9a237c9..2f3a27267 100644 --- a/python/cuspatial/cuspatial/__init__.py +++ b/python/cuspatial/cuspatial/__init__.py @@ -10,6 +10,7 @@ pairwise_point_distance, pairwise_point_linestring_distance, pairwise_point_linestring_nearest_points, + pairwise_point_polygon_distance, point_in_polygon, points_in_spatial_window, polygon_bounding_boxes, diff --git a/python/cuspatial/cuspatial/_lib/cpp/distance/point_polygon_distance.pxd b/python/cuspatial/cuspatial/_lib/cpp/distance/point_polygon_distance.pxd new file mode 100644 index 000000000..020bd4574 --- /dev/null +++ b/python/cuspatial/cuspatial/_lib/cpp/distance/point_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/point_polygon_distance.hpp" \ + namespace "cuspatial" nogil: + cdef unique_ptr[column] pairwise_point_polygon_distance( + const geometry_column_view & multipoints, + const geometry_column_view & multipolygons + ) except + diff --git a/python/cuspatial/cuspatial/_lib/distance.pyx b/python/cuspatial/cuspatial/_lib/distance.pyx index 9319278ee..f3d68717e 100644 --- a/python/cuspatial/cuspatial/_lib/distance.pyx +++ b/python/cuspatial/cuspatial/_lib/distance.pyx @@ -1,12 +1,15 @@ -# Copyright (c) 2022, NVIDIA CORPORATION. +# Copyright (c) 2022-2023, NVIDIA CORPORATION. -from libcpp.memory cimport unique_ptr +from libcpp.memory cimport make_shared, shared_ptr, unique_ptr from libcpp.utility cimport move from cudf._lib.column cimport Column from cudf._lib.cpp.column.column cimport column from cudf._lib.cpp.column.column_view cimport column_view +from cuspatial._lib.cpp.column.geometry_column_view cimport ( + geometry_column_view, +) from cuspatial._lib.cpp.distance.linestring_distance cimport ( pairwise_linestring_distance as c_pairwise_linestring_distance, ) @@ -16,7 +19,12 @@ from cuspatial._lib.cpp.distance.point_distance cimport ( from cuspatial._lib.cpp.distance.point_linestring_distance cimport ( pairwise_point_linestring_distance as c_pairwise_point_linestring_distance, ) +from cuspatial._lib.cpp.distance.point_polygon_distance cimport ( + pairwise_point_polygon_distance as c_pairwise_point_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 from cuspatial._lib.utils cimport unwrap_pyoptcol @@ -102,3 +110,35 @@ def pairwise_point_linestring_distance( c_linestring_points_xy, )) return Column.from_unique_ptr(move(c_result)) + + +def pairwise_point_polygon_distance( + point_collection_type, + Column multipoints, + Column multipolygons +): + + cdef collection_type_id point_multi_type = collection_type_py_to_c( + point_collection_type + ) + + cdef shared_ptr[geometry_column_view] c_multipoints = \ + make_shared[geometry_column_view]( + multipoints.view(), + point_multi_type, + geometry_type_id.POINT) + + 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_point_polygon_distance( + c_multipoints.get()[0], c_multipolygons.get()[0] + )) + + return Column.from_unique_ptr(move(c_result)) diff --git a/python/cuspatial/cuspatial/_lib/types.pxd b/python/cuspatial/cuspatial/_lib/types.pxd index 710724f47..6d3d1b3f9 100644 --- a/python/cuspatial/cuspatial/_lib/types.pxd +++ b/python/cuspatial/cuspatial/_lib/types.pxd @@ -2,6 +2,13 @@ from libc.stdint cimport uint8_t +from cuspatial._lib.cpp.types cimport collection_type_id, geometry_type_id + ctypedef uint8_t underlying_geometry_type_id_t ctypedef uint8_t underlying_collection_type_id_t + + +cdef geometry_type_id geometry_type_py_to_c(typ) except* + +cdef collection_type_id collection_type_py_to_c(typ) except* diff --git a/python/cuspatial/cuspatial/_lib/types.pyx b/python/cuspatial/cuspatial/_lib/types.pyx index 03354b5f1..103342cb2 100644 --- a/python/cuspatial/cuspatial/_lib/types.pyx +++ b/python/cuspatial/cuspatial/_lib/types.pyx @@ -28,3 +28,10 @@ class CollectionType(IntEnum): MULTI = ( collection_type_id.MULTI ) + + +cdef geometry_type_id geometry_type_py_to_c(typ : GeometryType): + return ( (typ.value)) + +cdef collection_type_id collection_type_py_to_c(typ : CollectionType): + return ( (typ.value)) diff --git a/python/cuspatial/cuspatial/core/spatial/__init__.py b/python/cuspatial/cuspatial/core/spatial/__init__.py index b269707d6..ee07838f9 100644 --- a/python/cuspatial/cuspatial/core/spatial/__init__.py +++ b/python/cuspatial/cuspatial/core/spatial/__init__.py @@ -1,4 +1,4 @@ -# Copyright (c) 2022, NVIDIA CORPORATION. +# Copyright (c) 2022-2023, NVIDIA CORPORATION. from .bounding import linestring_bounding_boxes, polygon_bounding_boxes from .distance import ( @@ -7,6 +7,7 @@ pairwise_linestring_distance, pairwise_point_distance, pairwise_point_linestring_distance, + pairwise_point_polygon_distance, ) from .filtering import points_in_spatial_window from .indexing import quadtree_on_points @@ -26,6 +27,7 @@ "sinusoidal_projection", "pairwise_point_distance", "pairwise_linestring_distance", + "pairwise_point_polygon_distance", "pairwise_point_linestring_distance", "pairwise_point_linestring_nearest_points", "polygon_bounding_boxes", diff --git a/python/cuspatial/cuspatial/core/spatial/distance.py b/python/cuspatial/cuspatial/core/spatial/distance.py index 9bf9afcdb..5d61e0564 100644 --- a/python/cuspatial/cuspatial/core/spatial/distance.py +++ b/python/cuspatial/cuspatial/core/spatial/distance.py @@ -10,16 +10,19 @@ pairwise_linestring_distance as cpp_pairwise_linestring_distance, 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, ) from cuspatial._lib.hausdorff import ( directed_hausdorff_distance as cpp_directed_hausdorff_distance, ) from cuspatial._lib.spatial import haversine_distance as cpp_haversine_distance +from cuspatial._lib.types import CollectionType from cuspatial.core.geoseries import GeoSeries from cuspatial.utils.column_utils import ( contains_only_linestrings, contains_only_multipoints, contains_only_points, + contains_only_polygons, ) @@ -382,6 +385,104 @@ def pairwise_point_linestring_distance( ) +def pairwise_point_polygon_distance(points: GeoSeries, polygons: GeoSeries): + """Compute distance between pairs of (multi)points and (multi)polygons + + The distance between a (multi)point and a (multi)polygon + is defined as the shortest distance between every point in the + multipoint and every edge of the (multi)polygon. If the multipoint 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)point in `points` and the ith + (multi)polygon in `polygons`. + + Parameters + ---------- + points : GeoSeries + The (multi)points 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, `points` series cannot contain both points and polygons. + Currently, it is unsupported that `points` contains both points and + multipoints. + + Examples + -------- + Compute distance between a point and a polygon: + >>> from shapely.geometry import Point + >>> points = cuspatial.GeoSeries([Point(0, 0)]) + >>> polygons = cuspatial.GeoSeries([Point(1, 1).buffer(0.5)]) + >>> cuspatial.pairwise_point_polygon_distance(points, polygons) + 0 0.914214 + dtype: float64 + + Compute distance between a multipoint and a multipolygon + + >>> from shapely.geometry import MultiPoint + >>> mpoints = cuspatial.GeoSeries([MultiPoint([Point(0, 0), Point(1, 1)])]) + >>> mpolys = cuspatial.GeoSeries([ + ... MultiPoint([Point(2, 2), Point(1, 2)]).buffer(0.5)]) + >>> cuspatial.pairwise_point_polygon_distance(mpoints, mpolys) + 0 0.5 + dtype: float64 + """ + + if len(points) != len(polygons): + raise ValueError("Unmatched input geoseries length.") + + if len(points) == 0: + return cudf.Series(dtype=points.points.xy.dtype) + + if not contains_only_points(points): + raise ValueError("`points` array must contain only points") + + if not contains_only_polygons(polygons): + raise ValueError("`linestrings` array must contain only linestrings") + + if len(points.points.xy) > 0 and len(points.multipoints.xy) > 0: + raise NotImplementedError( + "Mixing point and multipoint geometries is not supported" + ) + + point_collection_type = ( + CollectionType.SINGLE + if len(points.points.xy > 0) + else CollectionType.MULTI + ) + + # Handle slicing in geoseries + points_column = ( + points._column.points._column + if point_collection_type == CollectionType.SINGLE + else points._column.mpoints._column + ) + points_column = points_column.take( + points._column._meta.union_offsets._column + ) + + polygon_column = polygons._column.polygons._column + polygon_column = polygon_column.take( + polygons._column._meta.union_offsets._column + ) + + return Series._from_data( + { + None: c_pairwise_point_polygon_distance( + point_collection_type, points_column, polygon_column + ) + } + ) + + def _flatten_point_series( points: GeoSeries, ) -> Tuple[ diff --git a/python/cuspatial/cuspatial/tests/conftest.py b/python/cuspatial/cuspatial/tests/conftest.py index 4cb492e53..1c37cee77 100644 --- a/python/cuspatial/cuspatial/tests/conftest.py +++ b/python/cuspatial/cuspatial/tests/conftest.py @@ -1,4 +1,4 @@ -# Copyright (c) 2020-2021, NVIDIA CORPORATION. +# Copyright (c) 2020-2023, NVIDIA CORPORATION. import geopandas as gpd import numpy as np @@ -305,3 +305,13 @@ def factory(length): return mask return factory + + +@pytest.fixture +def naturalearth_cities(): + return gpd.read_file(gpd.datasets.get_path("naturalearth_cities")) + + +@pytest.fixture +def naturalearth_lowres(): + return gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) 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 new file mode 100644 index 000000000..199c41208 --- /dev/null +++ b/python/cuspatial/cuspatial/tests/spatial/distance/test_pairwise_point_polygon_distance.py @@ -0,0 +1,142 @@ +import geopandas as gpd +import pytest +from shapely.geometry import MultiPoint, MultiPolygon, Point, Polygon + +import cudf +from cudf.testing import assert_series_equal + +import cuspatial + + +def test_point_polygon_empty(): + lhs = cuspatial.GeoSeries.from_points_xy([]) + rhs = cuspatial.GeoSeries.from_polygons_xy([], [0], [0], [0]) + + got = cuspatial.pairwise_point_polygon_distance(lhs, rhs) + + expect = cudf.Series([], dtype="f8") + + assert_series_equal(got, expect) + + +def test_multipoint_polygon_empty(): + lhs = cuspatial.GeoSeries.from_multipoints_xy([], [0]) + rhs = cuspatial.GeoSeries.from_polygons_xy([], [0], [0], [0]) + + got = cuspatial.pairwise_point_polygon_distance(lhs, rhs) + + expect = cudf.Series([], dtype="f8") + + assert_series_equal(got, expect) + + +@pytest.mark.parametrize( + "points", [[Point(0, 0)], [MultiPoint([(1, 1), (2, 2)])]] +) +@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(points, polygons): + lhs = gpd.GeoSeries(points) + rhs = gpd.GeoSeries(polygons) + + dlhs = cuspatial.GeoSeries(points) + drhs = cuspatial.GeoSeries(polygons) + + expect = lhs.distance(rhs) + got = cuspatial.pairwise_point_polygon_distance(dlhs, drhs) + + assert_series_equal(got, cudf.Series(expect)) + + +@pytest.mark.parametrize( + "points", + [ + [Point(0, 0), Point(3, -3)], + [MultiPoint([(1, 1), (2, 2)]), MultiPoint([(3, 3), (4, 4)])], + ], +) +@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(points, polygons): + lhs = gpd.GeoSeries(points) + rhs = gpd.GeoSeries(polygons) + + dlhs = cuspatial.GeoSeries(points) + drhs = cuspatial.GeoSeries(polygons) + + expect = lhs.distance(rhs) + got = cuspatial.pairwise_point_polygon_distance(dlhs, drhs) + + assert_series_equal(got, cudf.Series(expect)) + + +def test_point_polygon_large(point_generator, polygon_generator): + N = 100 + points = gpd.GeoSeries(point_generator(N)) + polygons = gpd.GeoSeries(polygon_generator(N, 1.0, 1.5)) + + dpoints = cuspatial.from_geopandas(points) + dpolygons = cuspatial.from_geopandas(polygons) + + expect = points.distance(polygons) + got = cuspatial.pairwise_point_polygon_distance(dpoints, dpolygons) + + assert_series_equal(got, cudf.Series(expect)) + + +def test_point_polygon_geocities(naturalearth_cities, naturalearth_lowres): + N = 100 + 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] + ) + + got = cuspatial.pairwise_point_polygon_distance( + gpu_cities[:N], gpu_countries[:N] + ) + + assert_series_equal(cudf.Series(expect), got)