agri.field-area

export function fieldArea(boundary: readonly GeoPoint[]): FieldArea

/** The enclosed area on a sphere of radius 6,371,008.8 m, three ways. */
export interface FieldArea {
  /** to the whole square metre */
  readonly squareMetres: number;
  /** to 0.0001 ha (1 m²) */
  readonly hectares: number;
  /** international acres (4046.8564224 m²), to 0.0001 */
  readonly acres: number;
}

import { fieldArea } from "#fune/agri.field-area@^1";

import { type GeoPoint } from "./geo_point_in_polygon_types.ts";
import { sinCos } from "./math_sin_cos.ts";  ← from math.sin-cos ^1.0.0 · built alongside by fune
import { roundFloat } from "./math_round_float.ts";  ← from math.round-float ^1.0.0 · built alongside by fune
import { type FieldArea } from "./agri_field_area_types.ts";

const RADIUS = 6371008.8;
const RADIANS = 3.141592653589793 / 180;
const ACRE = 4046.8564224;

/**
 * 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.
 */
export function fieldArea(boundary: readonly GeoPoint[]): FieldArea {
  boundary.forEach((p, i) => {
    if (typeof p.lat !== "number" || !Number.isFinite(p.lat) || p.lat < -90 || p.lat > 90) {
      throw new RangeError(`point ${i + 1}: latitude must be a finite number from -90 to 90 degrees`);
    }
    if (typeof p.lng !== "number" || !Number.isFinite(p.lng) || p.lng < -180 || p.lng > 180) {
      throw new RangeError(`point ${i + 1}: longitude must be a finite number from -180 to 180 degrees`);
    }
  });
  let ring = boundary;
  if (ring.length > 1 && ring[0].lat === ring[ring.length - 1].lat && ring[0].lng === ring[ring.length - 1].lng) {
    ring = ring.slice(0, -1);
  }
  const distinct = new Set(ring.map((p) => `${p.lat},${p.lng}`));
  if (distinct.size < 3) throw new RangeError("a boundary needs at least 3 distinct points");

  const n = ring.length;
  const x: number[] = [];
  const y: number[] = [];
  const sin0 = sinCos(ring[0].lat * RADIANS).sin;
  let previous = 0;
  for (let i = 0; i < n; i++) {
    let d = ring[i].lng - 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.push(d * RADIANS);
    y.push(sinCos(ring[i].lat * RADIANS).sin - sin0);
  }
  let sum = 0;
  for (let i = 0; i < n; i++) {
    sum += (x[(i + 1) % n] - x[(i + n - 1) % n]) * y[i];
  }
  const area = (Math.abs(sum) * RADIUS * RADIUS) / 2;
  return {
    squareMetres: roundFloat(area, 0),
    hectares: roundFloat(area / 10000, 4),
    acres: roundFloat(area / ACRE, 4),
  };
}

Install

fune build

With that line in your source, in a TypeScript project (language typescript in fune.project), fune build resolves it and its 3 dependencies, pins them in fune.lock, downloads only the TypeScript 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 TypeScript implementation. Install it without the registry with fune add ./agri.field-area-1.0.0-typescript.fune, or fetch it from a terminal with fune pull agri.field-area@1.0.0:typescript.

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