Functional Weave
Code in Rust

agri.field-area@1.0.0

README.md

2,887 bytes · view raw

# agri.field-area

The area of a field from its boundary: the corners as latitude/longitude
points (`GeoPoint`, from `geo.point-in-polygon`), in order, walked either way.
The answer is in square metres, hectares and acres. A field near Cambridge
with corners 0.0009° of latitude and 0.0014° of longitude apart is 9548 m²,
0.9548 ha, 2.3595 acres.

## Method

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.