Add centroid, boundary and convex_hull function (returning new Geogra…

…phy) (#20)

Co-authored-by: Benoit Bovy <[email protected]>

CMakeLists.txt

@@ -69,6 +69,7 @@ endif()
 
 add_library(spherely MODULE
   src/geography.cpp
+  src/accessors-geog.cpp
   src/predicates.cpp
   src/spherely.cpp
   )

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

src/accessors-geog.cpp

@@ -0,0 +1,72 @@
+#include <s2geography.h>
+
+#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> point =
+        make_geography<s2geog::PointGeography, spherely::Point>(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<spherely::Geography>(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<spherely::Polygon>(std::move(s2_obj));
+    return PyObjectGeography::from_geog(std::move(geog_ptr));
+}
+
+void init_accessors(py::module& m) {
+    m.def("centroid",
+          py::vectorize(&centroid),
+          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");
+}

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 <class T1, class T2, class S>
-std::unique_ptr<T2> make_geography(S &&s2_obj) {
-    S2GeographyPtr s2geog_ptr = std::make_unique<T1>(std::forward<S>(s2_obj));
-    return std::make_unique<T2>(std::move(s2geog_ptr));
-}
-
 class PointFactory {
 public:
     static std::unique_ptr<Point> FromLatLonDegrees(double lat_degrees, double lon_degrees) {

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 <class T1, class T2, class S>
+std::unique_ptr<T2> make_geography(S&& s2_obj) {
+    S2GeographyPtr s2geog_ptr = std::make_unique<T1>(std::forward<S>(s2_obj));
+    return std::make_unique<T2>(std::move(s2geog_ptr));
+}
+
 }  // namespace spherely
 
 #endif  // SPHERELY_GEOGRAPHY_H_

src/pybind11.hpp

@@ -85,8 +85,10 @@ class PyObjectGeography : public py::object {
     // move semantics (Python takes ownership).
     //
     template <class T, std::enable_if_t<std::is_base_of<Geography, T>::value, bool> = true>
-    static py::object as_py_object(std::unique_ptr<T> geog_ptr) {
-        return py::cast(std::move(geog_ptr));
+    static PyObjectGeography from_geog(std::unique_ptr<T> geog_ptr) {
+        auto pyobj = py::cast(std::move(geog_ptr));
+        auto pyobj_geog = static_cast<PyObjectGeography &>(pyobj);
+        return std::move(pyobj_geog);
     }
 
     // Just check whether the object is a Geography

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);

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]: ...

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)