ENH: add area function

src/accessors-geog.cpp

@@ -32,6 +32,10 @@ double distance(PyObjectGeography a, PyObjectGeography b, double radius = EARTH_
     return s2geog::s2_distance(a_index, b_index) * radius;
 }
 
+double area(PyObjectGeography a, double radius = EARTH_RADIUS_METERS) {
+    return s2geog::s2_area(a.as_geog_ptr()->geog()) * radius * radius;
+}
+
 void init_accessors(py::module& m) {
     m.attr("EARTH_RADIUS_METERS") = py::float_(EARTH_RADIUS_METERS);
 
@@ -92,4 +96,20 @@ void init_accessors(py::module& m) {
             Radius of Earth in meters, default 6,371,010
 
     )pbdoc");
+
+    m.def("area",
+          py::vectorize(&area),
+          py::arg("a"),
+          py::arg("radius") = EARTH_RADIUS_METERS,
+          R"pbdoc(
+        Calculate the area of the geography.
+
+        Parameters
+        ----------
+        a : :py:class:`Geography` or array_like
+            Geography object
+        radius : float, optional
+            Radius of Earth in meters, default 6,371,010
+
+    )pbdoc");
 }

src/spherely.pyi

@@ -114,6 +114,18 @@ class _VFunc_Nin2optradius_Nout1(
         self, a: npt.ArrayLike, b: Geography, radius: float = ...
     ) -> npt.NDArray[_ArrayReturnDType]: ...
 
+class _VFunc_Nin1optradius_Nout1(
+    Generic[_NameType, _ScalarReturnType, _ArrayReturnDType]
+):
+    @property
+    def __name__(self) -> _NameType: ...
+    @overload
+    def __call__(self, a: Geography, radius: float = ...) -> _ScalarReturnType: ...
+    @overload
+    def __call__(
+        self, a: npt.ArrayLike, radius: float = ...
+    ) -> npt.NDArray[_ArrayReturnDType]: ...
+
 # Geography properties
 
 get_dimensions: _VFunc_Nin1_Nout1[Literal["get_dimensions"], Geography, Any]
@@ -190,6 +202,7 @@ convex_hull: _VFunc_Nin1_Nout1[
     Literal["convex_hull"], PolygonGeography, PolygonGeography
 ]
 distance: _VFunc_Nin2optradius_Nout1[Literal["distance"], float, float]
+area: _VFunc_Nin1optradius_Nout1[Literal["area"], float, float]
 
 # io functions
 

tests/test_accessors.py

@@ -1,3 +1,5 @@
+import math
+
 import numpy as np
 import pytest
 
@@ -124,3 +126,36 @@ def test_distance_with_custom_radius() -> None:
     )
     assert isinstance(actual, float)
     assert actual == pytest.approx(np.pi / 2)
+
+
+def test_area():
+    # scalar
+    geog = spherely.polygon([(0, 0), (90, 0), (0, 90), (0, 0)])
+    result = spherely.area(geog, radius=1)
+    assert isinstance(result, float)
+    expected = 4 * math.pi / 8
+    assert result == pytest.approx(expected, 1e-9)
+
+    result = spherely.area(geog)
+    assert result == pytest.approx(expected * spherely.EARTH_RADIUS_METERS**2, 1e-9)
+
+    # array
+    actual = spherely.area([geog], radius=1)
+    assert isinstance(actual, np.ndarray)
+    actual = actual[0]
+    assert isinstance(actual, float)
+    assert actual == pytest.approx(4 * math.pi / 8, 1e-9)
+
+
+@pytest.mark.parametrize(
+    "geog",
+    [
+        "POINT (-64 45)",
+        "POINT EMPTY",
+        "LINESTRING (0 0, 1 1)",
+        "LINESTRING EMPTY",
+        "POLYGON EMPTY",
+    ],
+)
+def test_area_empty(geog):
+    assert spherely.area(spherely.from_wkt(geog)) == 0