diff --git a/src/creation.cpp b/src/creation.cpp index 0ee2480..16a48b3 100644 --- a/src/creation.cpp +++ b/src/creation.cpp @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -64,26 +65,23 @@ std::vector make_s2points(const std::vector& points) { return std::move(s2points); } -// create a S2Loop from coordinates or Point objects. +void drop_closed_ring_duplicate(std::vector& s2points) { + if (!s2points.empty() && s2points.front() == s2points.back()) { + s2points.pop_back(); + } +} + +// create a S2Loop from a pre-built vector of S2Point. Callers are responsible +// for dropping any closed-ring duplicate vertex before calling. // // Normalization (to CCW order for identifying the loop interior) and validation // are both enabled by default. // -// Additional normalization is made here: -// - if the input loop is already closed, remove one of the end nodes -// // TODO: add option to skip normalization. // -template -std::unique_ptr make_s2loop(const std::vector& vertices, - bool check = true, - bool oriented = false) { - auto s2points = make_s2points(vertices); - - if (s2points.front() == s2points.back()) { - s2points.pop_back(); - } - +std::unique_ptr make_s2loop_from_points(const std::vector& s2points, + bool check = true, + bool oriented = false) { auto loop_ptr = std::make_unique(); loop_ptr->set_s2debug_override(S2Debug::DISABLE); @@ -101,13 +99,23 @@ std::unique_ptr make_s2loop(const std::vector& vertices, loop_ptr->Normalize(); } - return std::move(loop_ptr); + return loop_ptr; +} + +template +std::unique_ptr make_s2loop(const std::vector& vertices, + bool check = true, + bool oriented = false) { + auto s2points = make_s2points(vertices); + drop_closed_ring_duplicate(s2points); + return make_s2loop_from_points(s2points, check, oriented); } // create a S2Polygon. // std::unique_ptr make_s2polygon(std::vector> loops, - bool oriented = false) { + bool oriented = false, + bool check = true) { auto polygon_ptr = std::make_unique(); polygon_ptr->set_s2debug_override(S2Debug::DISABLE); @@ -118,7 +126,7 @@ std::unique_ptr make_s2polygon(std::vector> l } // Note: this also checks each loop of the polygon - if (!polygon_ptr->IsValid()) { + if (check && !polygon_ptr->IsValid()) { std::stringstream err; S2Error s2err; err << "polygon is not valid: "; @@ -159,6 +167,44 @@ py::array_t points(const py::array_t& coords) { return points; } +py::array_t polygons(const py::array_t& shells, + bool oriented, + bool validate) { + auto shells_data = shells.unchecked<3>(); + if (shells_data.shape(2) != 2) { + throw std::runtime_error("shells array must be of shape (N, K, 2)"); + } + + auto npolys = shells_data.shape(0); + auto nverts = shells_data.shape(1); + + auto result = py::array_t(npolys); + py::buffer_info buf = result.request(); + py::object* data = static_cast(buf.ptr); + + // Scratch buffer reused across polygons. Per-ring validity is always + // skipped; bad rings still surface via the polygon-level IsValid() in + // make_s2polygon when validate=true. + std::vector pts; + pts.reserve(static_cast(nverts)); + + for (py::ssize_t i = 0; i < npolys; i++) { + pts.clear(); + for (py::ssize_t j = 0; j < nverts; j++) { + pts.push_back(make_s2point(shells_data(i, j, 0), shells_data(i, j, 1))); + } + drop_closed_ring_duplicate(pts); + + std::vector> loops; + loops.push_back(make_s2loop_from_points(pts, /*check=*/false, oriented)); + + data[i] = make_py_geography( + make_s2polygon(std::move(loops), oriented, validate)); + } + + return result; +} + template std::unique_ptr create_multipoint(const std::vector& pts) { try { @@ -513,4 +559,51 @@ void init_creation(py::module& m) { A array of longitude, latitude coordinate tuples (i.e., with shape (N, 2)). )pbdoc"); + + m.def("polygons", + &polygons, + py::arg("shells"), + py::arg("oriented") = false, + py::arg("validate") = true, + R"pbdoc(polygons(shells, oriented=False, validate=True) + + Create an array of polygons from a numpy array of ring coordinates. + + This vectorized constructor is functionally equivalent to calling + :py:func:`create_polygon` for each shell but avoids the per-polygon + Python parsing overhead, making it much faster when building many + polygons of uniform shape (e.g. grid cells). + + Limitations compared to :py:func:`create_polygon`: all polygons + must have the same number of shell vertices (the numpy ndarray + input enforces this uniformity), and holes are not supported — + use :py:func:`create_polygon` in a Python loop if either varies + across the batch. + + Parameters + ---------- + shells : array_like + Array of shape ``(N, K, 2)`` giving ``N`` shell rings, each with + ``K`` vertices expressed as ``(longitude, latitude)`` in degrees. + Rings may be open (first vertex not repeated as last) or closed; + in the latter case the duplicate closing vertex is dropped. + 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). + By default (False), each ring is normalized so that the polygon + corresponds to the smaller area on the sphere. + validate : bool, default True + Validate each resulting polygon via ``S2Polygon::IsValid``. Set + to False only if the input polygons are known to be valid (e.g. + rectilinear grid cells); skipping the check is a meaningful + speedup on large batches but an invalid polygon will then be + constructed silently. + + Returns + ------- + polygons : ndarray + A 1-d object array of shape ``(N,)`` containing POLYGON + :class:`~spherely.Geography` objects. + + )pbdoc"); } diff --git a/src/spherely.pyi b/src/spherely.pyi index e9f03d0..dd8d9f2 100644 --- a/src/spherely.pyi +++ b/src/spherely.pyi @@ -212,6 +212,11 @@ def create_collection(geographies: Iterable[Geography]) -> GeometryCollection: . def points( longitude: npt.ArrayLike, latitude: npt.ArrayLike ) -> PointGeography | T_NDArray_Geography: ... +def polygons( + shells: npt.NDArray[np.float64], + oriented: bool = False, + validate: bool = True, +) -> T_NDArray_Geography: ... # Geography utils diff --git a/tests/test_creation.py b/tests/test_creation.py index 5c3c948..07b727f 100644 --- a/tests/test_creation.py +++ b/tests/test_creation.py @@ -1,3 +1,4 @@ +import numpy as np import pytest import spherely @@ -387,6 +388,97 @@ def test_multipolygon_invalid_geography() -> None: spherely.create_multipolygon([poly, line]) +def test_polygons_vectorized_basic() -> None: + shells = np.array( + [ + [(0, 0), (2, 0), (2, 2), (0, 2)], + [(4, 0), (6, 0), (6, 2), (4, 2)], + [(10, -5), (12, -5), (12, -3), (10, -3)], + [(-30, 40), (-28, 40), (-28, 42), (-30, 42)], + ], + dtype=np.float64, + ) + result = spherely.polygons(shells) + + assert isinstance(result, np.ndarray) + assert result.shape == (4,) + assert result.dtype == object + assert repr(result[0]).startswith("POLYGON ((0 0") + assert repr(result[1]).startswith("POLYGON ((4 0") + + # each vectorized polygon matches its scalar create_polygon equivalent + scalar = [spherely.create_polygon(ring.tolist()) for ring in shells] + assert [spherely.to_wkt(p) for p in result] == [spherely.to_wkt(p) for p in scalar] + + +def test_polygons_vectorized_closed_ring() -> None: + # open and closed input rings should produce the same polygon + open_ring = np.array( + [[(0, 0), (2, 0), (2, 2), (0, 2)]], + dtype=np.float64, + ) + closed_ring = np.array( + [[(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)]], + dtype=np.float64, + ) + a = spherely.polygons(open_ring)[0] + b = spherely.polygons(closed_ring)[0] + assert spherely.to_wkt(a) == spherely.to_wkt(b) + + +def test_polygons_vectorized_oriented() -> None: + # same pattern as test_polygon_oriented but using the vectorized path + shells = np.array( + [[(0, 0), (0, 2), (2, 2), (2, 0)]], + dtype=np.float64, + ) + (poly_cw,) = spherely.polygons(shells, oriented=True) + + point = spherely.create_point(1, 1) + # CW polygon + oriented=True => point at (1, 1) is NOT in the interior + assert not spherely.contains(poly_cw, point) + assert repr(poly_cw) == "POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))" + + +def test_polygons_vectorized_shape_errors() -> None: + # wrong trailing dimension + with pytest.raises(RuntimeError, match="shape"): + spherely.polygons(np.zeros((2, 4, 3), dtype=np.float64)) + + # wrong number of dimensions (surfaces from pybind11's unchecked view) + with pytest.raises(ValueError, match="number of dimensions"): + spherely.polygons(np.zeros((4, 2), dtype=np.float64)) + + +def test_polygons_vectorized_invalid_ring() -> None: + # same validation behavior as scalar create_polygon: bad rings surface as + # "polygon is not valid" from the final polygon-level check. + shells = np.array( + [ + [(0, 0), (2, 0), (2, 2), (0, 2)], + [(0, 0), (0, 2), (0, 2), (2, 0)], # invalid: duplicate + self-cross + ], + dtype=np.float64, + ) + with pytest.raises(ValueError, match="polygon is not valid"): + spherely.polygons(shells) + + +def test_polygons_vectorized_validate_false_skips_validation() -> None: + # same bad shell as above, but validate=False must not raise. + shells = np.array( + [ + [(0, 0), (2, 0), (2, 2), (0, 2)], + [(0, 0), (0, 2), (0, 2), (2, 0)], + ], + dtype=np.float64, + ) + result = spherely.polygons(shells, validate=False) + # we still get a Geography back for each input; no guarantee it's valid. + assert result.shape == (2,) + assert all(isinstance(g, spherely.Geography) for g in result) + + def test_collection() -> None: objs = [ spherely.create_point(0, 0),