ENH: add boolean operations (overlays) (#62)

CMakeLists.txt

@@ -73,6 +73,7 @@ endif()
 
 set(CPP_SOURCES
   src/accessors-geog.cpp
+  src/boolean-operations.cpp
   src/creation.cpp
   src/geography.cpp
   src/io.cpp

docs/api.rst

@@ -84,6 +84,19 @@ Predicates
    covers
    covered_by
 
+.. _api_overlays:
+
+Overlays (boolean operations)
+-----------------------------
+
+.. autosummary::
+   :toctree: _api_generated/
+
+   union
+   intersection
+   difference
+   symmetric_difference
+
 .. _api_constructive_ops:
 
 Constructive operations

src/boolean-operations.cpp

@@ -0,0 +1,91 @@
+#include <s2geography.h>
+#include <s2geography/geography.h>
+
+#include "constants.hpp"
+#include "creation.hpp"
+#include "geography.hpp"
+#include "pybind11.hpp"
+
+namespace py = pybind11;
+namespace s2geog = s2geography;
+using namespace spherely;
+
+class BooleanOp {
+public:
+    BooleanOp(S2BooleanOperation::OpType op_type)
+        : m_op_type(op_type), m_options(s2geog::GlobalOptions()) {
+        // TODO make this configurable
+        // m_options.polyline_layer_action = s2geography::GlobalOptions::OUTPUT_ACTION_IGNORE;
+    }
+
+    PyObjectGeography operator()(PyObjectGeography a, PyObjectGeography b) const {
+        const auto& a_index = a.as_geog_ptr()->geog_index();
+        const auto& b_index = b.as_geog_ptr()->geog_index();
+        std::unique_ptr<s2geog::Geography> geog_out =
+            s2geog::s2_boolean_operation(a_index, b_index, m_op_type, m_options);
+
+        return make_py_geography(std::move(geog_out));
+    }
+
+private:
+    S2BooleanOperation::OpType m_op_type;
+    s2geog::GlobalOptions m_options;
+};
+
+void init_boolean_operations(py::module& m) {
+    m.def("union",
+          py::vectorize(BooleanOp(S2BooleanOperation::OpType::UNION)),
+          py::arg("a"),
+          py::arg("b"),
+          R"pbdoc(
+        Computes the union of both geographies.
+
+        Parameters
+        ----------
+        a, b : :py:class:`Geography` or array_like
+            Geography object
+
+    )pbdoc");
+
+    m.def("intersection",
+          py::vectorize(BooleanOp(S2BooleanOperation::OpType::INTERSECTION)),
+          py::arg("a"),
+          py::arg("b"),
+          R"pbdoc(
+        Computes the intersection of both geographies.
+
+        Parameters
+        ----------
+        a, b : :py:class:`Geography` or array_like
+            Geography object
+
+    )pbdoc");
+
+    m.def("difference",
+          py::vectorize(BooleanOp(S2BooleanOperation::OpType::DIFFERENCE)),
+          py::arg("a"),
+          py::arg("b"),
+          R"pbdoc(
+        Computes the difference of both geographies.
+
+        Parameters
+        ----------
+        a, b : :py:class:`Geography` or array_like
+            Geography object
+
+    )pbdoc");
+
+    m.def("symmetric_difference",
+          py::vectorize(BooleanOp(S2BooleanOperation::OpType::SYMMETRIC_DIFFERENCE)),
+          py::arg("a"),
+          py::arg("b"),
+          R"pbdoc(
+        Computes the symmetric difference of both geographies.
+
+        Parameters
+        ----------
+        a, b : :py:class:`Geography` or array_like
+            Geography object
+
+    )pbdoc");
+}

src/spherely.cpp

@@ -8,6 +8,7 @@ namespace py = pybind11;
 void init_geography(py::module&);
 void init_creation(py::module&);
 void init_predicates(py::module&);
+void init_boolean_operations(py::module&);
 void init_accessors(py::module&);
 void init_io(py::module&);
 #if defined(S2GEOGRAPHY_VERSION_MAJOR) && \
@@ -28,6 +29,7 @@ PYBIND11_MODULE(spherely, m) {
     init_geography(m);
     init_creation(m);
     init_predicates(m);
+    init_boolean_operations(m);
     init_accessors(m);
     init_io(m);
 #if defined(S2GEOGRAPHY_VERSION_MAJOR) && \

src/spherely.pyi

@@ -211,6 +211,15 @@ touches: _VFunc_Nin2_Nout1[Literal["touches"], bool, bool]
 covers: _VFunc_Nin2_Nout1[Literal["covers"], bool, bool]
 covered_by: _VFunc_Nin2_Nout1[Literal["covered_by"], bool, bool]
 
+# boolean operations
+
+union: _VFunc_Nin2_Nout1[Literal["union"], Geography, Geography]
+intersection: _VFunc_Nin2_Nout1[Literal["intersection"], Geography, Geography]
+difference: _VFunc_Nin2_Nout1[Literal["difference"], Geography, Geography]
+symmetric_difference: _VFunc_Nin2_Nout1[
+    Literal["symmetric_difference"], Geography, Geography
+]
+
 # coords
 
 get_x: _VFunc_Nin1_Nout1[Literal["get_x"], float, np.float64]

tests/test_boolean_operations.py

@@ -0,0 +1,149 @@
+from packaging.version import Version
+
+import pytest
+
+import spherely
+
+poly1 = spherely.from_wkt("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))")
+poly2 = spherely.from_wkt("POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))")
+
+
+@pytest.mark.parametrize(
+    "geog1, geog2, expected",
+    [
+        ("POINT (30 10)", "POINT EMPTY", "POINT (30 10)"),
+        ("POINT EMPTY", "POINT EMPTY", "GEOMETRYCOLLECTION EMPTY"),
+        (
+            "LINESTRING (-45 0, 0 0)",
+            "LINESTRING (0 0, 0 10)",
+            "LINESTRING (-45 0, 0 0, 0 10)",
+        ),
+    ],
+)
+def test_union(geog1, geog2, expected):
+    result = spherely.union(spherely.from_wkt(geog1), spherely.from_wkt(geog2))
+    assert str(result) == expected
+
+
+def test_union_polygon():
+    result = spherely.union(poly1, poly2)
+
+    expected_near = (
+        spherely.area(spherely.from_wkt("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))"))
+        + spherely.area(spherely.from_wkt("POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))"))
+        - spherely.area(spherely.from_wkt("POLYGON ((5 5, 10 5, 10 15, 5 10, 5 5))"))
+    )
+    pytest.approx(spherely.area(result), expected_near)
+
+
+@pytest.mark.parametrize(
+    "geog1, geog2, expected",
+    [
+        ("POINT (30 10)", "POINT (30 10)", "POINT (30 10)"),
+        (
+            "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))",
+            "LINESTRING (0 5, 10 5)",
+            "LINESTRING (0 5, 10 5)",
+        ),
+    ],
+)
+def test_intersection(geog1, geog2, expected):
+    result = spherely.intersection(spherely.from_wkt(geog1), spherely.from_wkt(geog2))
+    assert str(result) == expected
+
+
+def test_intersection_empty():
+    result = spherely.intersection(poly1, spherely.from_wkt("POLYGON EMPTY"))
+    # assert spherely.is_empty(result)
+    assert str(result) == "GEOMETRYCOLLECTION EMPTY"
+
+    result = spherely.intersection(spherely.from_wkt("POLYGON EMPTY"), poly1)
+    assert str(result) == "GEOMETRYCOLLECTION EMPTY"
+
+    result = spherely.intersection(
+        spherely.from_wkt("POINT (0 1)"), spherely.from_wkt("POINT (1 2)")
+    )
+    assert str(result) == "GEOMETRYCOLLECTION EMPTY"
+
+
+def test_intersection_lines():
+    result = spherely.intersection(
+        spherely.from_wkt("LINESTRING (-45 0, 45 0)"),
+        spherely.from_wkt("LINESTRING (0 -10, 0 10)"),
+    )
+    assert str(result) == "POINT (0 0)"
+    assert spherely.distance(result, spherely.from_wkt("POINT (0 0)")) == 0
+
+
+def test_intersection_polygons():
+    result = spherely.intersection(poly1, poly2)
+    # TODO precision could be higher with snap level
+    precision = 2 if Version(spherely.__s2geography_version__) < Version("0.2.0") else 1
+    assert (
+        spherely.to_wkt(result, precision=precision)
+        == "POLYGON ((5 5, 10 5, 10 10, 5 10, 5 5))"
+    )
+
+
+def test_intersection_polygon_model():
+    poly = spherely.from_wkt("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))")
+    point = spherely.from_wkt("POINT (0 0)")
+
+    result = spherely.intersection(poly, point)
+    assert str(result) == "GEOMETRYCOLLECTION EMPTY"
+
+    # TODO this will be different depending on the model
+
+
+@pytest.mark.parametrize(
+    "geog1, geog2, expected",
+    [
+        ("POINT (30 10)", "POINT EMPTY", "POINT (30 10)"),
+        ("POINT EMPTY", "POINT EMPTY", "GEOMETRYCOLLECTION EMPTY"),
+        (
+            "LINESTRING (0 0, 45 0)",
+            "LINESTRING (0 0, 45 0)",
+            "GEOMETRYCOLLECTION EMPTY",
+        ),
+    ],
+)
+def test_difference(geog1, geog2, expected):
+    result = spherely.difference(spherely.from_wkt(geog1), spherely.from_wkt(geog2))
+    assert spherely.equals(result, spherely.from_wkt(expected))
+
+
+def test_difference_polygons():
+    result = spherely.difference(poly1, poly2)
+    expected_near = spherely.area(poly1) - spherely.area(
+        spherely.from_wkt("POLYGON ((5 5, 10 5, 10 10, 5 10, 5 5))")
+    )
+    pytest.approx(spherely.area(result), expected_near)
+
+
+@pytest.mark.parametrize(
+    "geog1, geog2, expected",
+    [
+        ("POINT (30 10)", "POINT EMPTY", "POINT (30 10)"),
+        ("POINT (30 10)", "POINT (30 10)", "GEOMETRYCOLLECTION EMPTY"),
+        ("POINT (30 10)", "POINT (30 20)", "MULTIPOINT ((30 20), (30 10))"),
+        (
+            "LINESTRING (0 0, 45 0)",
+            "LINESTRING (0 0, 45 0)",
+            "GEOMETRYCOLLECTION EMPTY",
+        ),
+    ],
+)
+def test_symmetric_difference(geog1, geog2, expected):
+    result = spherely.symmetric_difference(
+        spherely.from_wkt(geog1), spherely.from_wkt(geog2)
+    )
+    assert spherely.equals(result, spherely.from_wkt(expected))
+
+
+def test_symmetric_difference_polygons():
+    result = spherely.symmetric_difference(poly1, poly2)
+    expected_near = 2 * (
+        spherely.area(poly1)
+        - spherely.area(spherely.from_wkt("POLYGON ((5 5, 10 5, 10 10, 5 10, 5 5))"))
+    )
+    pytest.approx(spherely.area(result), expected_near)