From d40a8d25fe6d18bd65d713cf9d5c658668c1a728 Mon Sep 17 00:00:00 2001 From: Joris Van den Bossche Date: Fri, 29 Nov 2024 08:32:19 +0100 Subject: [PATCH] ENH: add boolean operations (overlays) (#62) --- CMakeLists.txt | 1 + docs/api.rst | 13 +++ src/boolean-operations.cpp | 91 +++++++++++++++++++ src/spherely.cpp | 2 + src/spherely.pyi | 9 ++ tests/test_boolean_operations.py | 149 +++++++++++++++++++++++++++++++ 6 files changed, 265 insertions(+) create mode 100644 src/boolean-operations.cpp create mode 100644 tests/test_boolean_operations.py diff --git a/CMakeLists.txt b/CMakeLists.txt index a9e0bb1..550b67f 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/docs/api.rst b/docs/api.rst index cc03070..42938bf 100644 --- a/docs/api.rst +++ b/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 diff --git a/src/boolean-operations.cpp b/src/boolean-operations.cpp new file mode 100644 index 0000000..aad4cb6 --- /dev/null +++ b/src/boolean-operations.cpp @@ -0,0 +1,91 @@ +#include +#include + +#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 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"); +} diff --git a/src/spherely.cpp b/src/spherely.cpp index bd0bb66..5c8e7e7 100644 --- a/src/spherely.cpp +++ b/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) && \ diff --git a/src/spherely.pyi b/src/spherely.pyi index 25dbcee..a663d6c 100644 --- a/src/spherely.pyi +++ b/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] diff --git a/tests/test_boolean_operations.py b/tests/test_boolean_operations.py new file mode 100644 index 0000000..8c4cc2c --- /dev/null +++ b/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)