impl/typescript.ts

Imports name this capability’s declared dependencies, which fune builds next to it in your project; each one links to its page.

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