Added more missing predicates

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

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
 

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