From 38d3bf9df81e86140dacb962b573fc95e44f9316 Mon Sep 17 00:00:00 2001 From: Joris Van den Bossche Date: Mon, 19 Aug 2024 16:24:29 +0200 Subject: [PATCH] Add centroid, boundary and convex_hull function (returning new Geography) (#20) Co-authored-by: Benoit Bovy --- CMakeLists.txt | 1 + docs/api.rst | 12 ++++++ src/accessors-geog.cpp | 72 +++++++++++++++++++++++++++++++++++ src/geography.cpp | 13 ------- src/geography.hpp | 13 +++++++ src/pybind11.hpp | 6 ++- src/spherely.cpp | 2 + src/spherely.pyi | 6 +++ tests/test_accessors.py | 83 +++++++++++++++++++++++++++++++++++++++++ 9 files changed, 193 insertions(+), 15 deletions(-) create mode 100644 src/accessors-geog.cpp create mode 100644 tests/test_accessors.py diff --git a/CMakeLists.txt b/CMakeLists.txt index cd5b3f4..fd9b2e9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -69,6 +69,7 @@ endif() add_library(spherely MODULE src/geography.cpp + src/accessors-geog.cpp src/predicates.cpp src/spherely.cpp ) diff --git a/docs/api.rst b/docs/api.rst index a8a9d2d..ea7ea46 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -42,6 +42,18 @@ Geography creation prepare destroy_prepared +.. _api_accessors: + +Geography accessors +------------------- + +.. autosummary:: + :toctree: _api_generated/ + + centroid + boundary + convex_hull + .. _api_predicates: Predicates diff --git a/src/accessors-geog.cpp b/src/accessors-geog.cpp new file mode 100644 index 0000000..172b786 --- /dev/null +++ b/src/accessors-geog.cpp @@ -0,0 +1,72 @@ +#include + +#include "geography.hpp" +#include "pybind11.hpp" + +namespace py = pybind11; +namespace s2geog = s2geography; +using namespace spherely; + +PyObjectGeography centroid(PyObjectGeography a) { + const auto& a_ptr = a.as_geog_ptr()->geog(); + auto s2_point = s2geog::s2_centroid(a_ptr); + std::unique_ptr point = + make_geography(s2_point); + return PyObjectGeography::from_geog(std::move(point)); +} + +PyObjectGeography boundary(PyObjectGeography a) { + const auto& a_ptr = a.as_geog_ptr()->geog(); + auto s2_obj = s2geog::s2_boundary(a_ptr); + // TODO return specific subclass + auto geog_ptr = std::make_unique(std::move(s2_obj)); + return PyObjectGeography::from_geog(std::move(geog_ptr)); +} + +PyObjectGeography convex_hull(PyObjectGeography a) { + const auto& a_ptr = a.as_geog_ptr()->geog(); + auto s2_obj = s2geog::s2_convex_hull(a_ptr); + auto geog_ptr = std::make_unique(std::move(s2_obj)); + return PyObjectGeography::from_geog(std::move(geog_ptr)); +} + +void init_accessors(py::module& m) { + m.def("centroid", + py::vectorize(¢roid), + py::arg("a"), + R"pbdoc( + Computes the centroid of each geography. + + Parameters + ---------- + a : :py:class:`Geography` or array_like + Geography object + + )pbdoc"); + + m.def("boundary", + py::vectorize(&boundary), + py::arg("a"), + R"pbdoc( + Computes the boundary of each geography. + + Parameters + ---------- + a : :py:class:`Geography` or array_like + Geography object + + )pbdoc"); + + m.def("convex_hull", + py::vectorize(&convex_hull), + py::arg("a"), + R"pbdoc( + Computes the convex hull of each geography. + + Parameters + ---------- + a : :py:class:`Geography` or array_like + Geography object + + )pbdoc"); +} diff --git a/src/geography.cpp b/src/geography.cpp index a165cde..1556092 100644 --- a/src/geography.cpp +++ b/src/geography.cpp @@ -46,19 +46,6 @@ S2Point to_s2point(const Point *vertex) { return vertex->s2point(); } -/* -** Helper to create Geography object wrappers. -** -** @tparam T1 The S2Geography wrapper type -** @tparam T2 This library wrapper type. -** @tparam S The S2Geometry type -*/ -template -std::unique_ptr make_geography(S &&s2_obj) { - S2GeographyPtr s2geog_ptr = std::make_unique(std::forward(s2_obj)); - return std::make_unique(std::move(s2geog_ptr)); -} - class PointFactory { public: static std::unique_ptr FromLatLonDegrees(double lat_degrees, double lon_degrees) { diff --git a/src/geography.hpp b/src/geography.hpp index 7826b3a..b1936d9 100644 --- a/src/geography.hpp +++ b/src/geography.hpp @@ -113,6 +113,19 @@ class Polygon : public Geography { } }; +/* +** Helper to create Geography object wrappers. +** +** @tparam T1 The S2Geography wrapper type +** @tparam T2 This library wrapper type. +** @tparam S The S2Geometry type +*/ +template +std::unique_ptr make_geography(S&& s2_obj) { + S2GeographyPtr s2geog_ptr = std::make_unique(std::forward(s2_obj)); + return std::make_unique(std::move(s2geog_ptr)); +} + } // namespace spherely #endif // SPHERELY_GEOGRAPHY_H_ diff --git a/src/pybind11.hpp b/src/pybind11.hpp index afd95bc..fc172dd 100644 --- a/src/pybind11.hpp +++ b/src/pybind11.hpp @@ -85,8 +85,10 @@ class PyObjectGeography : public py::object { // move semantics (Python takes ownership). // template ::value, bool> = true> - static py::object as_py_object(std::unique_ptr geog_ptr) { - return py::cast(std::move(geog_ptr)); + static PyObjectGeography from_geog(std::unique_ptr geog_ptr) { + auto pyobj = py::cast(std::move(geog_ptr)); + auto pyobj_geog = static_cast(pyobj); + return std::move(pyobj_geog); } // Just check whether the object is a Geography diff --git a/src/spherely.cpp b/src/spherely.cpp index 12c196e..df1d957 100644 --- a/src/spherely.cpp +++ b/src/spherely.cpp @@ -7,6 +7,7 @@ namespace py = pybind11; void init_geography(py::module&); void init_predicates(py::module&); +void init_accessors(py::module&); PYBIND11_MODULE(spherely, m) { m.doc() = R"pbdoc( @@ -19,6 +20,7 @@ PYBIND11_MODULE(spherely, m) { init_geography(m); init_predicates(m); + init_accessors(m); #ifdef VERSION_INFO m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); diff --git a/src/spherely.pyi b/src/spherely.pyi index 64ea736..425280c 100644 --- a/src/spherely.pyi +++ b/src/spherely.pyi @@ -107,6 +107,12 @@ contains: _VFunc_Nin2_Nout1[Literal["contains"], bool, bool] within: _VFunc_Nin2_Nout1[Literal["within"], bool, bool] disjoint: _VFunc_Nin2_Nout1[Literal["disjoint"], bool, bool] +# geography accessors + +centroid: _VFunc_Nin1_Nout1[Literal["centroid"], Geography, Point] +boundary: _VFunc_Nin1_Nout1[Literal["boundary"], Geography, Geography] +convex_hull: _VFunc_Nin1_Nout1[Literal["convex_hull"], Geography, Polygon] + # temp (remove) def create(arg0: Iterable[float], arg1: Iterable[float]) -> npt.NDArray[Any]: ... diff --git a/tests/test_accessors.py b/tests/test_accessors.py new file mode 100644 index 0000000..205823a --- /dev/null +++ b/tests/test_accessors.py @@ -0,0 +1,83 @@ +import numpy as np + +import spherely + +import pytest + + +@pytest.mark.parametrize( + "geog, expected", + [ + (spherely.Point(0, 0), spherely.Point(0, 0)), + ( + spherely.LineString([(0, 0), (0, 2)]), + spherely.Point(0, 1), + ), + ( + spherely.Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]), + spherely.Point(1, 1), + ), + ], +) +def test_centroid(geog, expected) -> None: + # scalar + actual = spherely.centroid(geog) + assert isinstance(actual, spherely.Point) + # TODO add some way of testing almost equality + # assert spherely.equals(actual, expected) + + # array + actual = spherely.centroid([geog]) + assert isinstance(actual, np.ndarray) + actual = actual[0] + assert isinstance(actual, spherely.Point) + # assert spherely.equals(actual, expected) + + +@pytest.mark.parametrize( + "geog, expected", + [ + (spherely.Point(0, 0), "GEOMETRYCOLLECTION EMPTY"), + (spherely.LineString([(0, 0), (0, 2), (2, 2)]), "MULTIPOINT ((0 0), (2 2))"), + ( + spherely.Polygon([(0, 0), (2, 0), (2, 2), (1.5, 0.5)]), + "LINESTRING (0.5 1.5, 2 2, 0 2, 0 0, 0.5 1.5)", + ), + ], +) +def test_boundary(geog, expected) -> None: + # scalar + actual = spherely.boundary(geog) + assert str(actual) == expected + + # array + actual = spherely.boundary([geog]) + assert isinstance(actual, np.ndarray) + assert str(actual[0]) == expected + + +@pytest.mark.parametrize( + "geog, expected", + [ + ( + spherely.LineString([(0, 0), (0, 2), (2, 2)]), + spherely.Polygon([(0, 0), (0, 2), (2, 2)]), + ), + ( + spherely.Polygon([(0, 0), (2, 0), (2, 2), (1.5, 0.5)]), + spherely.Polygon([(0, 0), (2, 0), (2, 2)]), + ), + ], +) +def test_convex_hull(geog, expected) -> None: + # scalar + actual = spherely.convex_hull(geog) + assert isinstance(actual, spherely.Polygon) + assert spherely.equals(actual, expected) + + # array + actual = spherely.convex_hull([geog]) + assert isinstance(actual, np.ndarray) + actual = actual[0] + assert isinstance(actual, spherely.Polygon) + assert spherely.equals(actual, expected)