agri.field-area

def field_area(boundary: Sequence[GeoPoint]) -> FieldArea

@dataclass(frozen=True)
class FieldArea:
    """The enclosed area on a sphere of radius 6,371,008.8 m, three ways."""

    #: to the whole square metre
    square_metres: float
    #: to 0.0001 ha (1 m²)
    hectares: float
    #: international acres (4046.8564224 m²), to 0.0001
    acres: float

from fune.agri.field_area import field_area  # agri.field-area@^1

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),
    )

Install

fune build

With that line in your source, in a Python project (language python in fune.project), fune build resolves it and its 3 dependencies, pins them in fune.lock, downloads only the Python package of each, and builds the code above into your project’s .fune/build, one readable file per capability with a header linking back here. Or pin a range in fune.project and build in one step:

fune add agri.field-area

The manifest, vectors and README with only the Python implementation. Install it without the registry with fune add ./agri.field-area-1.0.0-python.fune, or fetch it from a terminal with fune pull agri.field-area@1.0.0:python.

The whole function, every language, is one file too: agri.field-area-1.0.0.fune, 20,067 bytes, sha256 2e2dc272023c9eff0d5dfb894cadecee64c19cb251b22efe036999e3cac3a1e7. It installs into a project of any language.

# fune: before agri.field-area

# fune: after agri.field-area

# fune: replace geo.point-in-polygon in agri.field-area
# fune: replace math.round-float in agri.field-area
# fune: replace math.sin-cos in agri.field-area

# fune: step agri.field-area after <n|label>

More from the author

The Earth is a sphere of radius 6,371,008.8 m, the IUGG mean radius that `geo.distance` also uses. The boundary is projected with the Lambert cylindrical equal-area projection (x = R x longitude in radians, y = R x sin latitude), which preserves area, and the enclosed area is the shoelace formula there:

A = R² / 2 x | sum over i of (lng[i+1] - lng[i-1]) x sin(lat[i]) |

This is the spherical polygon area formula of Chamberlain and Duquette (2007), used by many web mapping libraries. The edges are straight lines in that projection rather than great circles; for anything the size of a field the difference is far below the uncertainty of the corner positions.

**Accuracy.** A sphere is not the Earth: against the WGS 84 ellipsoid a spherical area is typically out by a few tenths of a percent (it depends on latitude), so a 10 ha field can be a few hundredths of a hectare out. That is fine for planning inputs and checking a sketch; claims for payment (Rural Payments, IACS) must use the official land parcel area, and survey work needs an ellipsoidal method.

## Determinism

Sines come from `math.sin-cos` (the platform maths libraries differ in the last bit between languages). Longitudes are taken relative to the first point and sines relative to the first point's sine, which keeps precision for small fields far from the origin, and every sum is formed in the same order in all three languages, so the unrounded area is the same double. Each result is then rounded half away from zero with `math.round-float`: square metres to 0 places, hectares and acres to 4 (acres are international acres of 4046.8564224 m²).

## Edge cases

- Closing the ring (repeating the first point last) is optional. - Direction does not matter; the absolute value is taken. - A field crossing the 180° meridian is measured the short way round: successive longitudes are unwrapped to within 180° of each other. - Collinear points enclose 0. - A self-intersecting boundary (a figure of eight) is not detected; its "area" is the difference of its loops. Holes are not supported: subtract them yourself. - Fewer than 3 distinct points, or a coordinate outside ±90 / ±180 (usually swapped latitude and longitude), is an error.

## Source

R. G. Chamberlain and W. H. Duquette, "Some Algorithms for Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion Laboratory (2007), https://trs.jpl.nasa.gov/handle/2014/41271. Earth radius: H. Moritz, "Geodetic Reference System 1980", Journal of Geodesy 74 (2000) 128-133.

Files