From ff7adee393aba2f6b1c6d01fe7374e76bd7081da Mon Sep 17 00:00:00 2001 From: Joel Jaeschke Date: Sun, 20 Oct 2024 17:19:12 +0200 Subject: [PATCH] Added more missing predicates --- src/predicates.cpp | 73 ++++++++++++++++++++++++++++++++++++++++ src/spherely.pyi | 3 ++ tests/test_predicates.py | 69 +++++++++++++++++++++++++++++++++++++ 3 files changed, 145 insertions(+) diff --git a/src/predicates.cpp b/src/predicates.cpp index 933b43a..fd5135d 100644 --- a/src/predicates.cpp +++ b/src/predicates.cpp @@ -117,4 +117,77 @@ void init_predicates(py::module& m) { Geography object(s) )pbdoc"); + + m.def("touches", + py::vectorize(Predicate([](const s2geog::ShapeIndexGeography& a_index, + const s2geog::ShapeIndexGeography& b_index, + const S2BooleanOperation::Options& options) { + S2BooleanOperation::Options closedOptions = options; + closedOptions.set_polygon_model(S2BooleanOperation::PolygonModel::CLOSED); + closedOptions.set_polyline_model(S2BooleanOperation::PolylineModel::CLOSED); + + S2BooleanOperation::Options openOptions = options; + openOptions.set_polygon_model(S2BooleanOperation::PolygonModel::OPEN); + openOptions.set_polyline_model(S2BooleanOperation::PolylineModel::OPEN); + + return s2geog::s2_intersects(a_index, b_index, closedOptions) && + !s2geog::s2_intersects(a_index, b_index, openOptions); + })), + py::arg("a"), + py::arg("b"), + R"pbdoc( + Returns True if A and B intersect, but their interiors do not intersect. + A and B have at least one point in common, where the common point + lies in at least one boundary. + + Parameters + ---------- + a, b : :py:class:`Geography` or array_like + Geography object(s) + + )pbdoc"); + + m.def("covers", + py::vectorize(Predicate([](const s2geog::ShapeIndexGeography& a_index, + const s2geog::ShapeIndexGeography& b_index, + const S2BooleanOperation::Options& options) { + S2BooleanOperation::Options closedOptions = options; + closedOptions.set_polyline_model(S2BooleanOperation::PolylineModel::CLOSED); + closedOptions.set_polygon_model(S2BooleanOperation::PolygonModel::CLOSED); + + return s2geog::s2_contains(a_index, b_index, closedOptions); + })), + py::arg("a"), + py::arg("b"), + R"pbdoc( + Returns True if every point in B lies inside the interior or boundary of A. + + Parameters + ---------- + a, b : :py:class:`Geography` or array_like + Geography object(s) + + )pbdoc"); + + m.def("covered_by", + py::vectorize(Predicate([](const s2geog::ShapeIndexGeography& a_index, + const s2geog::ShapeIndexGeography& b_index, + const S2BooleanOperation::Options& options) { + S2BooleanOperation::Options closedOptions = options; + closedOptions.set_polyline_model(S2BooleanOperation::PolylineModel::CLOSED); + closedOptions.set_polygon_model(S2BooleanOperation::PolygonModel::CLOSED); + + return s2geog::s2_contains(b_index, a_index, closedOptions); + })), + py::arg("a"), + py::arg("b"), + R"pbdoc( + Returns True if every point in A lies inside the interior or boundary of B. + + Parameters + ---------- + a, b : :py:class:`Geography` or array_like + Geography object(s) + + )pbdoc"); } diff --git a/src/spherely.pyi b/src/spherely.pyi index f01e74f..48be786 100644 --- a/src/spherely.pyi +++ b/src/spherely.pyi @@ -130,6 +130,9 @@ equals: _VFunc_Nin2_Nout1[Literal["intersects"], bool, bool] contains: _VFunc_Nin2_Nout1[Literal["contains"], bool, bool] within: _VFunc_Nin2_Nout1[Literal["within"], bool, bool] disjoint: _VFunc_Nin2_Nout1[Literal["disjoint"], bool, bool] +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] # geography accessors diff --git a/tests/test_predicates.py b/tests/test_predicates.py index ecf40d4..001dfea 100644 --- a/tests/test_predicates.py +++ b/tests/test_predicates.py @@ -100,3 +100,72 @@ def test_disjoint(): a2 = spherely.Point(50, 9) b2 = spherely.LineString([(50, 8), (60, 8)]) assert spherely.disjoint(a2, b2) + + +def test_touches(): + a = spherely.Polygon([(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)]) + b = np.array( + [ + spherely.Polygon([(1.0, 1.0), (1.0, 2.0), (2.0, 2.0), (2.0, 1.0)]), + spherely.Polygon([(0.5, 0.5), (0.5, 1.5), (1.5, 1.5), (1.5, 0.5)]), + ] + ) + + actual = spherely.touches(a, b) + expected = np.array([True, False]) + np.testing.assert_array_equal(actual, expected) + + a_p = spherely.Point(1.0, 1.0) + b_p = spherely.Point(1.0, 1.0) + # Points do not have a boundary, so they cannot touch per definition + # This is consistent with PostGIS for example (cmp. https://postgis.net/docs/ST_Touches.html) + assert not spherely.touches(a_p, b_p) + + b_line = spherely.LineString([(1.0, 1.0), (1.0, 2.0)]) + assert spherely.touches(a_p, b_line) + + +def test_covers(): + a = spherely.Polygon([(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)]) + b_p = np.array( + [ + spherely.Point(2.0, 2.0), + spherely.Point(1.0, 1.0), + spherely.Point(0.5, 0.5), + ] + ) + + actual = spherely.covers(a, b_p) + expected = np.array([False, True, True]) + np.testing.assert_array_equal(actual, expected) + + b_poly_in = spherely.Polygon([(0.1, 0.1), (0.1, 0.9), (0.9, 0.9), (0.9, 0.1)]) + b_poly_part_boundary = spherely.Polygon( + [(0.0, 0.0), (0.0, 0.75), (0.75, 0.75), (0.75, 0.0)] + ) + + assert spherely.covers(a, b_poly_in) + assert spherely.covers(a, b_poly_part_boundary) # This fails, but should not + + +def test_covered_by(): + a = spherely.Polygon([(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)]) + b_p = np.array( + [ + spherely.Point(2.0, 2.0), + spherely.Point(1.0, 1.0), + spherely.Point(0.5, 0.5), + ] + ) + + actual = spherely.covered_by(b_p, a) + expected = np.array([False, True, True]) + np.testing.assert_array_equal(actual, expected) + + b_poly_in = spherely.Polygon([(0.1, 0.1), (0.1, 0.9), (0.9, 0.9), (0.9, 0.1)]) + b_poly_part_boundary = spherely.Polygon( + [(0.0, 0.0), (0.0, 0.75), (0.75, 0.75), (0.75, 0.0)] + ) + + assert spherely.covered_by(b_poly_in, a) + assert spherely.covers(b_poly_part_boundary, a) # This fails, but should not