feat: implement shoelace formula for polygon area calculation

maths/bearing.py

@@ -0,0 +1,78 @@
+"""
+Initial bearing (forward azimuth) between two geographic points.
+
+Points are (latitude, longitude) in decimal degrees. Returns bearing in degrees
+clockwise from true north in the range [0, 360).
+
+Reference: https://en.wikipedia.org/wiki/Bearing_(navigation)
+"""
+
+from __future__ import annotations
+
+import math
+
+
+def initial_bearing(
+    origin: tuple[float, float], destination: tuple[float, float]
+) -> float:
+    """
+    Compute the initial bearing from `origin` to `destination`.
+
+    Parameters
+    ----------
+    origin, destination : tuple[float, float]
+        (latitude, longitude) in decimal degrees.
+
+    Returns
+    -------
+    float
+        Initial bearing in degrees, clockwise from north in [0, 360).
+
+    Raises
+    ------
+    TypeError
+        If inputs are not 2-tuples of numbers.
+    ValueError
+        If the two points are identical (bearing undefined).
+
+    Examples
+    >>> round(initial_bearing((50.066389, -5.714722), (58.643889, -3.07)), 3)
+    9.12
+    >>> round(initial_bearing((0.0, 0.0), (1.0, 1.0)), 3)
+    44.996
+    >>> initial_bearing((0.0, 0.0), (0.0, 0.0))
+    Traceback (most recent call last):
+    ...
+    ValueError: origin and destination are the same point; bearing is undefined
+    """
+    try:
+        lat1, lon1 = float(origin[0]), float(origin[1])
+        lat2, lon2 = float(destination[0]), float(destination[1])
+    except Exception as exc:
+        raise TypeError(
+            "origin and destination must be 2-tuples of numeric values"
+        ) from exc
+
+    if lat1 == lat2 and lon1 == lon2:
+        raise ValueError(
+            "origin and destination are the same point; bearing is undefined"
+        )
+
+    # convert degrees to radians
+    phi1 = math.radians(lat1)
+    phi2 = math.radians(lat2)
+    delta_lambda = math.radians(lon2 - lon1)
+
+    x = math.sin(delta_lambda) * math.cos(phi2)
+    y = math.cos(phi1) * math.sin(phi2) - math.sin(phi1) * math.cos(phi2) * math.cos(
+        delta_lambda
+    )
+
+    theta = math.atan2(x, y)  # result in radians relative to north
+    bearing = (math.degrees(theta) + 360.0) % 360.0
+    return bearing
+
+
+if __name__ == "__main__":
+    # simple demonstration
+    print(initial_bearing((50.066389, -5.714722), (58.643889, -3.07)))

maths/shoelace_area.py

@@ -0,0 +1,70 @@
+"""
+Shoelace formula (Gauss's area formula) for polygon area.
+
+The function accepts an iterable of (x, y) pairs and returns the polygon area
+as a non-negative float.
+
+References:
+- https://en.wikipedia.org/wiki/Shoelace_formula
+"""
+
+from __future__ import annotations
+
+from collections.abc import Iterable, Sequence
+
+
+def shoelace_area(points: Iterable[tuple[float, float]]) -> float:
+    """
+    Compute the area of a simple polygon using the shoelace formula.
+
+    Parameters
+    ----------
+    points:
+        Iterable of (x, y) coordinate pairs. Points may be ints or floats.
+        The polygon is assumed closed (the function will wrap the last point
+        to the first).
+
+    Returns
+    -------
+    float
+        Non-negative area of the polygon.
+
+    Raises
+    ------
+    ValueError
+        If fewer than 3 points are provided.
+    TypeError
+        If points are not pairs of numbers.
+
+    Examples
+    >>> shoelace_area([(0, 0), (4, 0), (0, 3)])
+    6.0
+    >>> shoelace_area([(0, 0), (1, 0), (1, 1), (0, 1)])
+    1.0
+    >>> shoelace_area(list(reversed([(0, 0), (2, 0), (2, 2), (0, 2)])))
+    4.0
+    >>> shoelace_area([(0, 0), (2, 0), (2, 2), (0, 2)])
+    4.0
+    """
+    pts = list(points)
+    n = len(pts)
+    if n < 3:
+        raise ValueError("At least 3 points are required to form a polygon")
+
+    try:
+        coords: Sequence[tuple[float, float]] = [(float(x), float(y)) for x, y in pts]
+    except Exception as exc:
+        raise TypeError("points must be an iterable of (x, y) numeric pairs") from exc
+
+    s = 0.0
+    for i in range(n):
+        x1, y1 = coords[i]
+        x2, y2 = coords[(i + 1) % n]
+        s += x1 * y2 - x2 * y1
+
+    return abs(s) / 2.0
+
+
+if __name__ == "__main__":
+    example = [(0, 0), (4, 0), (0, 3)]
+    print("example area:", shoelace_area(example))