feat: use precise and fast area intersection algorithm

package.json

@@ -44,6 +44,7 @@
     "test:integration": "vitest run --config vitest.integration.config.mts",
     "test:e2e": "vitest run --config vitest.e2e.config.mts",
     "test:coverage": "vitest run --coverage",
+    "bench": "vitest bench --run",
     "verify:fast": "npm run typecheck && npm run typecheck:test && npm run build && npm run test:unit",
     "verify": "npm run verify:fast && npm run test:integration",
     "verify:full": "npm run verify && npm run test:e2e",
@@ -56,16 +57,15 @@
   "dependencies": {
     "@ignfab/gpf-schema-store": "^0.2.2",
     "@rgrove/parse-xml": "^4.2.3",
-    "@turf/area": "^7.4.0",
     "@turf/bbox": "^7.3.5",
     "@turf/bbox-clip": "^7.4.0",
     "@turf/bbox-polygon": "^7.4.0",
     "@turf/centroid": "^7.3.5",
     "@turf/circle": "^7.4.0",
     "@turf/distance": "^7.3.5",
     "@turf/helpers": "^7.3.5",
-    "@turf/intersect": "^7.4.0",
     "@turf/length": "^7.4.0",
+    "earcut": "^3.2.4",
     "jsts": "^2.12.1",
     "mcp-framework": "^0.2.22",
     "node-fetch": "^3.3.2",
@@ -81,6 +81,7 @@
     "@langchain/ollama": "^1.3.0",
     "@modelcontextprotocol/inspector": "^2.0.0",
     "@modelcontextprotocol/sdk": "^1.30.0",
+    "@turf/intersect": "^7.4.0",
     "@types/geojson": "^7946.0.16",
     "@types/node": "^26.1.2",
     "@types/supertest": "^7.2.1",

src/helpers/area.ts

@@ -0,0 +1,50 @@
+import type { Geometry, Position } from "geojson";
+
+const DEGREES_TO_RADIANS = Math.PI / 180;
+export const EARTH_RADIUS = 6371008.7714;
+
+/** Area of a ring on the unit sphere, exact for edges straight in the lon/lat plane (GeoJSON standard lines).
+ *
+ * Unlike `@turf/area`, which approximates each edge, splitting an edge leaves it
+ * unchanged and zero-width spikes contribute nothing.
+ */
+export function sphericalRingArea(ring: Position[]) : number {
+  let sum = 0;
+  ring.forEach(([lon1, lat1], index) => {
+    const [lon2, lat2] = ring[(index + 1) % ring.length];
+    // Integral of sin(latitude) along the edge, over the longitude.
+    const halfDeltaLat = (lat2 - lat1) * DEGREES_TO_RADIANS / 2;
+    const sinc = halfDeltaLat == 0 ? 1 : Math.sin(halfDeltaLat) / halfDeltaLat;
+    sum += (lon2 - lon1) * DEGREES_TO_RADIANS * Math.sin((lat1 + lat2) * DEGREES_TO_RADIANS / 2) * sinc;
+  });
+  return Math.abs(sum);
+}
+
+/** Area of a geometry.
+ * 
+ * Replacement of @turf/area with improved precision.
+ * 
+ * @param geo Input geometry
+ * @returns The area of the areal parts of the geometry, in m²
+ */
+export default function area(geo: Geometry) : number {
+  switch (geo.type) {
+    case "Point":
+    case "MultiPoint":
+    case "LineString":
+    case "MultiLineString":
+      return 0;
+    case "Polygon":
+    case "MultiPolygon": {
+      let total = 0;
+      for (const poly of geo.type === "Polygon" ? [geo.coordinates] : geo.coordinates) {
+        const [outer, ...holes] = poly;
+        total += holes.reduce((remaining, hole) => remaining - sphericalRingArea(hole), sphericalRingArea(outer));
+      }
+      return total * EARTH_RADIUS**2;
+    }
+    case "GeometryCollection": {
+      return geo.geometries.reduce((acc, geom) => acc + area(geom), 0);
+    }
+  }
+}

src/helpers/geojson.ts

@@ -1,4 +1,4 @@
-import type { Geometry, GeometryCollection } from "geojson";
+import type { Geometry, GeometryCollection, MultiPolygon, Polygon } from "geojson";
 
 /**
  * Narrow guard for the minimal GeoJSON geometry shape `geometryToEwkt` requires.
@@ -15,3 +15,22 @@ export function isGeometryLike(value: unknown): value is Exclude<Geometry, Geome
     "coordinates" in value
   );
 }
+
+/** Strip the empty rings `@turf/bbox-clip` emits for parts lying outside the box.
+ *
+ * A fully-clipped-away Polygon comes back as `{ coordinates: [] }` and a
+ * MultiPolygon keeps one empty entry per discarded part (`[[...], []]`), both of
+ * which are invalid GeoJSON. Returns null when nothing areal survives.
+ */
+export function dropEmptyRings(geom: Geometry) : Polygon | MultiPolygon | null {
+  if (geom.type == "Polygon") {
+    return geom.coordinates.length == 0 ? null : geom;
+  }
+  if (geom.type == "MultiPolygon") {
+    const coordinates = geom.coordinates.filter((polygon) => polygon.length > 0);
+    return coordinates.length == 0 ? null : { type: "MultiPolygon", coordinates };
+  }
+  // `bboxClip` is typed over every geometry it accepts, but the caller only ever
+  // passes polygons, so a non-areal result means nothing areal survived.
+  return null;
+}