ENH: add oriented argument to create_polygon() (#85)

* create_polygon(): add oriented argument

* add multipolygon test

src/creation.cpp

@@ -78,7 +78,9 @@ std::vector<S2Point> make_s2points(const std::vector<V> &points) {
 // TODO: add option to skip normalization.
 //
 template <class V>
-std::unique_ptr<S2Loop> make_s2loop(const std::vector<V> &vertices, bool check = true) {
+std::unique_ptr<S2Loop> make_s2loop(const std::vector<V> &vertices,
+                                    bool check = true,
+                                    bool oriented = false) {
     auto s2points = make_s2points(vertices);
 
     if (s2points.front() == s2points.back()) {
@@ -98,17 +100,25 @@ std::unique_ptr<S2Loop> make_s2loop(const std::vector<V> &vertices, bool check =
         throw py::value_error(err.str());
     }
 
-    loop_ptr->Normalize();
+    if (!oriented) {
+        loop_ptr->Normalize();
+    }
 
     return std::move(loop_ptr);
 }
 
 // create a S2Polygon.
 //
-std::unique_ptr<S2Polygon> make_s2polygon(std::vector<std::unique_ptr<S2Loop>> loops) {
+std::unique_ptr<S2Polygon> make_s2polygon(std::vector<std::unique_ptr<S2Loop>> loops,
+                                          bool oriented = false) {
     auto polygon_ptr = std::make_unique<S2Polygon>();
     polygon_ptr->set_s2debug_override(S2Debug::DISABLE);
-    polygon_ptr->InitNested(std::move(loops));
+
+    if (oriented) {
+        polygon_ptr->InitOriented(std::move(loops));
+    } else {
+        polygon_ptr->InitNested(std::move(loops));
+    }
 
     // Note: this also checks each loop of the polygon
     if (!polygon_ptr->IsValid()) {
@@ -222,7 +232,8 @@ std::unique_ptr<Geography> create_multilinestring(const std::vector<Geography *>
 
 template <class V>
 std::unique_ptr<Geography> create_polygon(const std::vector<V> &shell,
-                                          const std::optional<std::vector<std::vector<V>>> &holes) {
+                                          const std::optional<std::vector<std::vector<V>>> &holes,
+                                          bool oriented = false) {
     // fastpath empty polygon
     if (shell.empty()) {
         if (holes.has_value() && !holes.value().empty()) {
@@ -234,18 +245,18 @@ std::unique_ptr<Geography> create_polygon(const std::vector<V> &shell,
     std::vector<std::unique_ptr<S2Loop>> loops;
 
     try {
-        loops.push_back(make_s2loop(shell, false));
+        loops.push_back(make_s2loop(shell, false, oriented));
     } catch (const EmptyGeographyException &error) {
         throw py::value_error("can't create Polygon with empty component");
     }
 
     if (holes.has_value()) {
         for (const auto &ring : holes.value()) {
-            loops.push_back(make_s2loop(ring, false));
+            loops.push_back(make_s2loop(ring, false, oriented));
         }
     }
 
-    return make_geography<s2geog::PolygonGeography>(make_s2polygon(std::move(loops)));
+    return make_geography<s2geog::PolygonGeography>(make_s2polygon(std::move(loops), oriented));
 }
 
 std::unique_ptr<Geography> create_multipolygon(const std::vector<Geography *> &polygons) {
@@ -369,15 +380,12 @@ void init_creation(py::module &m) {
 
     m.def(
          "create_polygon",
-         [](py::none, py::none) {
-             // TODO: remove explicit creation of S2Polygon, see
-             // https://github.com/paleolimbot/s2geography/pull/31
-             auto empty_poly = std::make_unique<S2Polygon>();
-             return make_geography(
-                 std::make_unique<s2geog::PolygonGeography>(std::move(empty_poly)));
+         [](py::none, py::none, bool) {
+             return make_geography(std::make_unique<s2geog::PolygonGeography>());
          },
          py::arg("shell") = py::none(),
          py::arg("holes") = py::none(),
+         py::arg("oriented") = false,
          R"pbdoc(create_polygon(shell: Sequence | None = None, holes: Sequence | None = None) -> Geography
         Create a POLYGON geography.
 
@@ -389,16 +397,24 @@ void init_creation(py::module &m) {
         holes : sequence, optional
             A list of sequences of objects where each sequence satisfies the same
             requirements as the ``shell`` argument.
+        oriented : bool, default False
+            Set to True if polygon ring directions are known to be correct
+            (i.e., shell ring vertices are defined counter clockwise and hole
+            ring vertices are defined clockwise).
+            By default (False), it will return the polygon with the smaller
+            area.
 
     )pbdoc")
         .def("create_polygon",
              &create_polygon<std::pair<double, double>>,
              py::arg("shell"),
-             py::arg("holes") = py::none())
+             py::arg("holes") = py::none(),
+             py::arg("oriented") = false)
         .def("create_polygon",
              &create_polygon<Geography *>,
              py::arg("shell"),
-             py::arg("holes") = py::none());
+             py::arg("holes") = py::none(),
+             py::arg("oriented") = false);
 
     m.def("create_multipolygon",
           &create_multipolygon,

src/spherely.pyi

@@ -190,16 +190,19 @@ def create_multilinestring(
 def create_polygon(
     shell: None = None,
     holes: None = None,
+    oriented: bool = False,
 ) -> PolygonGeography: ...
 @overload
 def create_polygon(
     shell: Iterable[Sequence[float]],
     holes: Iterable[Iterable[Sequence[float]]] | None = None,
+    oriented: bool = False,
 ) -> PolygonGeography: ...
 @overload
 def create_polygon(
     shell: Iterable[PointGeography],
     holes: Iterable[Iterable[PointGeography]] | None = None,
+    oriented: bool = False,
 ) -> PolygonGeography: ...
 def create_multipolygon(
     polygons: Iterable[PolygonGeography],

tests/test_creation.py

@@ -209,6 +209,33 @@ def test_polygon_normalize() -> None:
     assert repr(poly_cw) == "POLYGON ((2 0, 2 2, 0 2, 0 0, 2 0))"
 
 
+@pytest.mark.parametrize(
+    "shell",
+    [
+        [(0, 0), (0, 2), (2, 2), (2, 0)],
+        [
+            spherely.create_point(0, 0),
+            spherely.create_point(0, 2),
+            spherely.create_point(2, 2),
+            spherely.create_point(2, 0),
+        ],
+    ],
+)
+def test_polygon_oriented(shell) -> None:
+    # "CW" polygon + oriented=True => the polygon's interior is the largest
+    # area of the sphere divided by the polygon's ring
+    poly_cw = spherely.create_polygon(shell, oriented=True)
+
+    point = spherely.create_point(1, 1)
+
+    # point above is NOT in polygon's interior
+    assert not spherely.contains(poly_cw, point)
+
+    # polygon vertices not reordered
+    # TODO: better to test actual coordinate values when implemented
+    assert repr(poly_cw) == "POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))"
+
+
 @pytest.mark.parametrize(
     "shell,holes",
     [
@@ -271,6 +298,16 @@ def test_polygon_normalize_holes() -> None:
     )
 
 
+def test_polygon_oriented_holes() -> None:
+    # CCW polygon hole vertices + oriented=True => error
+    with pytest.raises(ValueError, match="Inconsistent loop orientations detected"):
+        spherely.create_polygon(
+            shell=[(0, 0), (2, 0), (2, 2), (0, 2)],
+            holes=[[(0.5, 0.5), (1.5, 0.5), (1.5, 1.5), (0.5, 1.5)]],
+            oriented=True,
+        )
+
+
 def test_polygon_empty() -> None:
     poly = spherely.create_polygon()
     assert repr(poly).startswith("POLYGON EMPTY")
@@ -340,6 +377,19 @@ def test_multipolygons() -> None:
     assert repr(multipoly).startswith("MULTIPOLYGON (((0 0")
 
 
+def test_multipolygons_oriented() -> None:
+    # same than `test_polygon_oriented`: make sure that it works for multipolygon too
+    poly_cw = spherely.create_polygon([(0, 0), (0, 2), (2, 2), (2, 0)], oriented=True)
+
+    # original polygon loops are cloned (so shouldn't be normalized) before being passed
+    # to the multipolygon
+    multipoly = spherely.create_multipolygon([poly_cw])
+
+    point = spherely.create_point(1, 1)
+
+    assert not spherely.contains(multipoly, point)
+
+
 def test_multipolygon_invalid_geography() -> None:
     poly = spherely.create_polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
     line = spherely.create_linestring([(3, 0), (3, 1)])