Imports name this capability’s declared dependencies, which fune builds next to it in your project; each one links to its page.
import math
from typing import List, Sequence
from .agri_field_area_types import FieldArea
from .geo_point_in_polygon_types import GeoPoint
from .math_round_float import round_float ← from math.round-float ^1.0.0 · built alongside by fune
from .math_sin_cos import sin_cos ← from math.sin-cos ^1.0.0 · built alongside by fune
RADIUS = 6371008.8
RADIANS = 3.141592653589793 / 180
ACRE = 4046.8564224
def _number(value: object) -> bool:
return isinstance(value, (int, float)) and not isinstance(value, bool) and math.isfinite(value)
def field_area(boundary: Sequence[GeoPoint]) -> FieldArea:
"""Area enclosed by a boundary, by the shoelace formula in the Lambert
cylindrical equal-area projection (x = R lng, y = R sin lat), which is the
spherical polygon formula of Chamberlain and Duquette (2007).
Longitudes are unwrapped from the first point and sines are taken relative
to the first point's, which keeps a one-hectare field accurate when the
coordinates are large; sin comes from math.sin-cos so every language adds
the same doubles in the same order.
"""
for i, p in enumerate(boundary):
if not _number(p.lat) or p.lat < -90 or p.lat > 90:
raise ValueError("point %d: latitude must be a finite number from -90 to 90 degrees" % (i + 1,))
if not _number(p.lng) or p.lng < -180 or p.lng > 180:
raise ValueError("point %d: longitude must be a finite number from -180 to 180 degrees" % (i + 1,))
ring = list(boundary)
if len(ring) > 1 and ring[0].lat == ring[-1].lat and ring[0].lng == ring[-1].lng:
ring = ring[:-1]
if len({(float(p.lat), float(p.lng)) for p in ring}) < 3:
raise ValueError("a boundary needs at least 3 distinct points")
n = len(ring)
x: List[float] = []
y: List[float] = []
sin0 = sin_cos(float(ring[0].lat) * RADIANS).sin
previous = 0.0
for p in ring:
d = float(p.lng) - float(ring[0].lng)
# Take the short way round, so a field across the antimeridian works.
while d - previous > 180:
d -= 360
while d - previous < -180:
d += 360
previous = d
x.append(d * RADIANS)
y.append(sin_cos(float(p.lat) * RADIANS).sin - sin0)
total = 0.0
for i in range(n):
total += (x[(i + 1) % n] - x[(i + n - 1) % n]) * y[i]
area = abs(total) * RADIUS * RADIUS / 2
return FieldArea(
square_metres=round_float(area, 0),
hectares=round_float(area / 10000, 4),
acres=round_float(area / ACRE, 4),
)