From 8e2d48105150f714040a54edfedf978da50d2865 Mon Sep 17 00:00:00 2001
From: Joel Jaeschke <joel.jaeschke@gmail.com>
Date: Tue, 26 Nov 2024 15:19:59 +0100
Subject: [PATCH] ENH: Added length and perimeter implementation (#77)

* Added length and (broken) perimeter implementation

* Added radius parameter to functions

* Add new functions to documentation
---
 docs/api.rst            |  2 ++
 src/accessors-geog.cpp  | 39 +++++++++++++++++++++++++++++++++
 src/spherely.pyi        |  2 ++
 tests/test_accessors.py | 48 +++++++++++++++++++++++++++++++++++++++++
 4 files changed, 91 insertions(+)

diff --git a/docs/api.rst b/docs/api.rst
index 08f7f3e..efa6696 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -62,6 +62,8 @@ Measurement
 
    area
    distance
+   length
+   perimeter
 
 .. _api_predicates:
 
diff --git a/src/accessors-geog.cpp b/src/accessors-geog.cpp
index 4d413df..27b7d6b 100644
--- a/src/accessors-geog.cpp
+++ b/src/accessors-geog.cpp
@@ -36,6 +36,14 @@ double area(PyObjectGeography a, double radius = EARTH_RADIUS_METERS) {
     return s2geog::s2_area(a.as_geog_ptr()->geog()) * radius * radius;
 }
 
+double length(PyObjectGeography a, double radius = EARTH_RADIUS_METERS) {
+    return s2geog::s2_length(a.as_geog_ptr()->geog()) * radius;
+}
+
+double perimeter(PyObjectGeography a, double radius = EARTH_RADIUS_METERS) {
+    return s2geog::s2_perimeter(a.as_geog_ptr()->geog()) * radius;
+}
+
 void init_accessors(py::module& m) {
     m.attr("EARTH_RADIUS_METERS") = py::float_(EARTH_RADIUS_METERS);
 
@@ -112,4 +120,35 @@ void init_accessors(py::module& m) {
             Radius of Earth in meters, default 6,371,010
 
     )pbdoc");
+
+    m.def("length",
+          py::vectorize(&length),
+          py::arg("a"),
+          py::arg("radius") = EARTH_RADIUS_METERS,
+          R"pbdoc(
+        Calculates the length of a line geography, returning zero for other types.
+
+        Parameters
+        ----------
+        a : :py:class:`Geography` or array_like
+            Geography object
+        radius : float, optional
+            Radius of Earth in meters, default 6,371,010
+
+   )pbdoc");
+
+    m.def("perimeter",
+          py::vectorize(&perimeter),
+          py::arg("a"),
+          py::arg("radius") = EARTH_RADIUS_METERS,
+          R"pbdoc(
+        Calculates the perimeter of a polygon geography, returning zero for other types.
+
+        Parameters
+        ----------
+        a : :py:class:`Geography` or array_like
+            Geography object
+        radius : float, optional
+            Radius of Earth in meters, default 6,371,010
+    )pbdoc");
 }
diff --git a/src/spherely.pyi b/src/spherely.pyi
index ea90227..b894673 100644
--- a/src/spherely.pyi
+++ b/src/spherely.pyi
@@ -213,6 +213,8 @@ convex_hull: _VFunc_Nin1_Nout1[
 ]
 distance: _VFunc_Nin2optradius_Nout1[Literal["distance"], float, float]
 area: _VFunc_Nin1optradius_Nout1[Literal["area"], float, float]
+length: _VFunc_Nin1optradius_Nout1[Literal["length"], float, float]
+perimeter: _VFunc_Nin1optradius_Nout1[Literal["perimeter"], float, float]
 
 # io functions
 
diff --git a/tests/test_accessors.py b/tests/test_accessors.py
index a203db3..319e8bc 100644
--- a/tests/test_accessors.py
+++ b/tests/test_accessors.py
@@ -159,3 +159,51 @@ def test_area():
 )
 def test_area_empty(geog):
     assert spherely.area(spherely.from_wkt(geog)) == 0
+
+
+def test_length():
+    geog = spherely.linestring([(0, 0), (1, 0)])
+    result = spherely.length(geog, radius=1)
+    assert isinstance(result, float)
+    expected = 1.0 * np.pi / 180.0
+    assert result == pytest.approx(expected, 1e-9)
+
+    actual = spherely.length([geog], radius=1)
+    assert isinstance(actual, np.ndarray)
+    actual = actual[0]
+    assert isinstance(actual, float)
+    assert actual == pytest.approx(expected, 1e-9)
+
+
+@pytest.mark.parametrize(
+    "geog",
+    [
+        "POINT (0 0)",
+        "POINT EMPTY",
+        "POLYGON EMPTY",
+        "POLYGON ((0 0, 0 1, 1 0, 0 0))",
+    ],
+)
+def test_length_invalid(geog):
+    assert spherely.length(spherely.from_wkt(geog)) == 0.0
+
+
+def test_perimeter():
+    geog = spherely.polygon([(0, 0), (0, 90), (90, 90), (90, 0), (0, 0)])
+    result = spherely.perimeter(geog, radius=1)
+    assert isinstance(result, float)
+    expected = 3 * 90 * np.pi / 180.0
+    assert result == pytest.approx(expected, 1e-9)
+
+    actual = spherely.perimeter([geog], radius=1)
+    assert isinstance(actual, np.ndarray)
+    actual = actual[0]
+    assert isinstance(actual, float)
+    assert actual == pytest.approx(expected, 1e-9)
+
+
+@pytest.mark.parametrize(
+    "geog", ["POINT (0 0)", "POINT EMPTY", "LINESTRING (0 0, 1 0)", "POLYGON EMPTY"]
+)
+def test_perimeter_invalid(geog):
+    assert spherely.perimeter(spherely.from_wkt(geog)) == 0.0