From 640878051711b346823501e39a93f2fc204a392e Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Mon, 28 Sep 2026 14:24:10 +0200 Subject: [PATCH 1/6] feat: use precise and fast area intersection algorithm --- package-lock.json | 81 ++++---------- package.json | 2 +- src/helpers/area.ts | 106 ++++++++++++++++++ src/wfs/spatialExtras.ts | 50 +++++++-- test/helpers/area.test.ts | 17 +++ test/wfs/response.test.ts | 225 +++++++++++++++++++++++++++++++++----- 6 files changed, 387 insertions(+), 94 deletions(-) create mode 100644 src/helpers/area.ts create mode 100644 test/helpers/area.test.ts diff --git a/package-lock.json b/package-lock.json index fb99e33b..a4c25ccf 100644 --- a/package-lock.json +++ b/package-lock.json @@ -11,7 +11,6 @@ "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", @@ -21,6 +20,7 @@ "@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", @@ -1620,48 +1620,6 @@ "dev": true, "license": "MIT" }, - "node_modules/@turf/area": { - "version": "7.4.0", - "resolved": "https://registry.npmjs.org/@turf/area/-/area-7.4.0.tgz", - "integrity": "sha512-B7q5f6QwKIxxbl4/L59lAYA+FF7h2hv6st96PuKue0TGCYGSxAG/lV8LHSPUSOuOrLywN4jytJcKK3vOFBs7Nw==", - "license": "MIT", - "dependencies": { - "@turf/helpers": "7.4.0", - "@turf/meta": "7.4.0", - "@types/geojson": "^7946.0.10", - "tslib": "^2.8.1" - }, - "funding": { - "url": "https://opencollective.com/turf" - } - }, - "node_modules/@turf/area/node_modules/@turf/helpers": { - "version": "7.4.0", - "resolved": "https://registry.npmjs.org/@turf/helpers/-/helpers-7.4.0.tgz", - "integrity": "sha512-7PAwLZqOdRzTI5g9bHvUwlloAXPDH/mlajtryk0tw4ZwGMtmXsAyF4QundsAMfy4u48Fyd3AUqMXY0SMxqJGWg==", - "license": "MIT", - "dependencies": { - "@types/geojson": "^7946.0.10", - "tslib": "^2.8.1" - }, - "funding": { - "url": "https://opencollective.com/turf" - } - }, - "node_modules/@turf/area/node_modules/@turf/meta": { - "version": "7.4.0", - "resolved": "https://registry.npmjs.org/@turf/meta/-/meta-7.4.0.tgz", - "integrity": "sha512-3cLUvlEyDuSnMSzrjhaLAEiYR8xhbfWyTVQlrlEw40xL81d4KF4PqUWbjTXKpXZStdYbet2GCurl79KMypNw6g==", - "license": "MIT", - "dependencies": { - "@turf/helpers": "7.4.0", - "@types/geojson": "^7946.0.10", - "tslib": "^2.8.1" - }, - "funding": { - "url": "https://opencollective.com/turf" - } - }, "node_modules/@turf/bbox": { "version": "7.3.5", "resolved": "https://registry.npmjs.org/@turf/bbox/-/bbox-7.3.5.tgz", @@ -3322,6 +3280,12 @@ "node": ">= 0.4" } }, + "node_modules/earcut": { + "version": "3.2.4", + "resolved": "https://registry.npmjs.org/earcut/-/earcut-3.2.4.tgz", + "integrity": "sha512-FHc6yNpnTORB0tM5kmZ03iD2mDuNakDE2lvf4Hw1XkfxtEHIGVi01B/GqUSwx+NcC1WGjidAfpTXzeALtPtfpg==", + "license": "ISC" + }, "node_modules/ecdsa-sig-formatter": { "version": "1.0.11", "resolved": "https://registry.npmjs.org/ecdsa-sig-formatter/-/ecdsa-sig-formatter-1.0.11.tgz", @@ -3626,9 +3590,9 @@ "license": "Unlicense" }, "node_modules/fast-uri": { - "version": "3.1.5", - "resolved": "https://registry.npmjs.org/fast-uri/-/fast-uri-3.1.5.tgz", - "integrity": "sha512-gHwA1O9LDIcKunMKhObS/HimwtehO1nPUECKAu5TpKgaO19fcWEl4bliWe1jWxVFvIXztJjjQ4L8XQ1EU9f7Jw==", + "version": "3.1.8", + "resolved": "https://registry.npmjs.org/fast-uri/-/fast-uri-3.1.8.tgz", + "integrity": "sha512-GZMtZUTNRpOVIECoXwLNZS5xUGE+mVNbTB8h/7Rwh2TFWcBQiPzTgyZi05BF9UMZKkLJv8XBRJTlU7zg8+ZfMg==", "funding": [ { "type": "github", @@ -3994,9 +3958,9 @@ } }, "node_modules/hono": { - "version": "4.13.1", - "resolved": "https://registry.npmjs.org/hono/-/hono-4.13.1.tgz", - "integrity": "sha512-kdJoFVv2xmayw6cY09H7AbMJMt8Jn5jdlEdXsP7AGBdF2DIptVlKlOLKXP41yPip4/a3yQPv9gVcJYI8YY04dw==", + "version": "4.13.10", + "resolved": "https://registry.npmjs.org/hono/-/hono-4.13.10.tgz", + "integrity": "sha512-dQuLsa5oO+47QVMVMaaD9cIv8ctmVtK1iRvwWngkfloFJMeFeuoUFDswIqZGxmGX3hrRzREgEArjkK9OgsQEhA==", "license": "MIT", "engines": { "node": ">=16.9.0" @@ -5676,12 +5640,13 @@ } }, "node_modules/qs": { - "version": "6.15.2", - "resolved": "https://registry.npmjs.org/qs/-/qs-6.15.2.tgz", - "integrity": "sha512-Rzq0KEyX/w/tEybncDgdkZrJgVUsUMk3xjh3t5bv3S1HTAtg+uOYt72+ZfwiQwKdysThkTBdL/rTi6HDmX9Ddw==", + "version": "6.16.0", + "resolved": "https://registry.npmjs.org/qs/-/qs-6.16.0.tgz", + "integrity": "sha512-h6fhOIaRrID2CbEY2fqs+7t+UXZo+MLAnU5gRIq85uFtdiUPCdsApMlHhXogKVM4HM2DVbIjGNTTYH2OcmP1vA==", "license": "BSD-3-Clause", "dependencies": { - "side-channel": "^1.1.0" + "es-define-property": "^1.0.1", + "side-channel": "^1.1.1" }, "engines": { "node": ">=0.6" @@ -5975,14 +5940,14 @@ } }, "node_modules/side-channel": { - "version": "1.1.0", - "resolved": "https://registry.npmjs.org/side-channel/-/side-channel-1.1.0.tgz", - "integrity": "sha512-ZX99e6tRweoUXqR+VBrslhda51Nh5MTQwou5tnUDgbtyM0dBgmhEDtWGP/xbKn6hqfPRHujUNwz5fy/wbbhnpw==", + "version": "1.1.1", + "resolved": "https://registry.npmjs.org/side-channel/-/side-channel-1.1.1.tgz", + "integrity": "sha512-6x6dK6zJdpTzF4sQeNYxwtvBzf6Eg4GtlesS94HOvTudUeyK2WXAaIfmDgsyslYrRBeFIlsi54AYsFGUuhmvrQ==", "license": "MIT", "dependencies": { "es-errors": "^1.3.0", - "object-inspect": "^1.13.3", - "side-channel-list": "^1.0.0", + "object-inspect": "^1.13.4", + "side-channel-list": "^1.0.1", "side-channel-map": "^1.0.1", "side-channel-weakmap": "^1.0.2" }, diff --git a/package.json b/package.json index 8f54f72e..f0450efb 100644 --- a/package.json +++ b/package.json @@ -56,7 +56,6 @@ "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", @@ -66,6 +65,7 @@ "@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", diff --git a/src/helpers/area.ts b/src/helpers/area.ts new file mode 100644 index 00000000..56ee1dad --- /dev/null +++ b/src/helpers/area.ts @@ -0,0 +1,106 @@ +import earcut, { deviation, flatten } from "earcut"; +import type { Geometry, Position } from "geojson"; + +type Triangle = [Position, Position, Position]; + +const DEGREES_TO_RADIANS = Math.PI / 180; +const EARTH_RADIUS = 6371008.7714; + +/** Positive when `c` lies on the left of the line going from `a` to `b`, negative on its right. */ +function cross(a: Position, b: Position, c: Position) : number { + return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]); +} + +/** Clip a closed ring by a triangle, using the Sutherland–Hodgman algorithm. + * + * The result may contain zero-width spikes along the triangle edges, which do not + * change its area. It is empty when nothing areal survives. + */ +export function clipRingToTriangle(ring: Position[], [a, b, c]: Triangle) : Position[] { + // Orient the triangle counterclockwise, so that its interior lies on the left of each edge. + const edges = cross(a, b, c) > 0 ? [[a, b], [b, c], [c, a]] : [[a, c], [c, b], [b, a]]; + let clipped = ring.slice(0, -1); // Drop the closing position + for (const [start, end] of edges) { + const input = clipped; + clipped = []; + input.forEach((current, index) => { + const previous = input.at(index - 1)!; // The last position for the first one + const previousSide = cross(start, end, previous); + const currentSide = cross(start, end, current); + if ((previousSide >= 0) !== (currentSide >= 0)) { + // The edge crosses the line: keep the crossing point. Both sides have + // different signs, so the denominator cannot be 0. + const t = previousSide / (previousSide - currentSide); + clipped.push([previous[0] + t * (current[0] - previous[0]), previous[1] + t * (current[1] - previous[1])]); + } + if (currentSide >= 0) { + clipped.push(current); + } + }); + } + return clipped.length < 3 ? [] : [...clipped, clipped[0]]; +} + +/** Split a polygon into triangles, or return null if it cannot be done reliably, + * as for self-intersecting rings. + */ +export function triangulate(polygon: Position[][]) : Triangle[] | null { + const { vertices, holes, dimensions } = flatten(polygon); + const indices = earcut(vertices, holes, dimensions); + // The triangles must cover the polygon exactly. + if (deviation(vertices, holes, dimensions, indices) > 1e-6) { + return null; + } + const position = (index: number) : Position => [vertices[index * dimensions], vertices[index * dimensions + 1]]; + const triangles : Triangle[] = []; + for (let i = 0; i < indices.length; i += 3) { + triangles.push([position(indices[i]), position(indices[i + 1]), position(indices[i + 2])]); + } + return triangles; +} + +/** 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); + } + } +} diff --git a/src/wfs/spatialExtras.ts b/src/wfs/spatialExtras.ts index 7a488f52..a490fc1a 100644 --- a/src/wfs/spatialExtras.ts +++ b/src/wfs/spatialExtras.ts @@ -1,7 +1,6 @@ import { centroid } from "@turf/centroid"; import { bbox } from "@turf/bbox"; import turfLength from "@turf/length"; -import { area } from "@turf/area"; import { intersect } from "@turf/intersect"; import { circle } from "@turf/circle"; import { bboxPolygon } from "@turf/bbox-polygon"; @@ -15,6 +14,7 @@ import type { SpatialFilter, } from "./schema.js"; import { bboxClip } from "@turf/bbox-clip"; +import area, { clipRingToTriangle, triangulate } from "../helpers/area.js"; export type FeatureCollectionPostProcessInput = { typename: string, @@ -190,6 +190,46 @@ export function dropEmptyRings(geom: Geometry) : Polygon | MultiPolygon | null { return null; } + +/** Return the area of the intersection between two areal (2D) geometries. + * + * Clipping by a triangle is simple, fast and robust, unlike general polygon + * clipping: `geo` is split into triangles, by which `filterPolygons` is clipped. + */ +function polygonsIntersectionArea(geo: Polygon | MultiPolygon, filterPolygons: Polygon | MultiPolygon) : number { + // Whatever lies outside the feature's bbox cannot intersect it, so clipping the + // reference first leaves the result unchanged while only the neighbouring + // vertices are processed afterwards. + const clippedFilter = dropEmptyRings(bboxClip(filterPolygons, bbox(geo)).geometry); + if (!clippedFilter) return 0; // no areal overlap + const filterParts = clippedFilter.type == "Polygon" ? [clippedFilter.coordinates] : clippedFilter.coordinates; + + let total = 0; + for (const polygon of geo.type == "Polygon" ? [geo.coordinates] : geo.coordinates) { + const polygonGeometry : Polygon = { type: "Polygon", coordinates: polygon }; + const triangles = triangulate(polygon); + if (!triangles) { + // Fall back on general polygon clipping. + const inter = intersect(featureCollection([feature(clippedFilter), feature(polygonGeometry)])); + total += inter == null ? 0 : area(inter.geometry); + continue; + } + let covered = 0; + for (const triangle of triangles) { + for (const [outer, ...holes] of filterParts) { + covered += area({ type: "Polygon", coordinates: [clipRingToTriangle(outer, triangle)]}); + for (const hole of holes) { + covered -= area({ type: "Polygon", coordinates: [clipRingToTriangle(hole, triangle)]}); + } + } + } + // Snap to 0 or to the whole polygon despite rounding errors. + const polygonArea = area(polygonGeometry); + total += covered > polygonArea * (1 - 1e-9) ? polygonArea : covered > polygonArea * 1e-9 ? covered : 0; + } + return total; +} + /** Return the area (m²) of the part of a geometry lying inside a spatial filter. * * null when it cannot be computed: the geometry or the filter has no areal part, @@ -210,13 +250,7 @@ function intersectionAreaWithSpatialFilter(geom: Geometry, spatialFilter: Spatia case "intersects_feature": case "travel_time": { if (!filterPolygons) return null; // non-areal filter, or filter preparation failed - // Whatever lies outside the feature's bbox cannot intersect it, so clipping the - // reference first leaves the result unchanged while polyclip only ever processes - // the neighbouring vertices. - const clippedFilter = dropEmptyRings(bboxClip(filterPolygons, bbox(geo)).geometry); - if (!clippedFilter) return 0; // no overlap - const inter = intersect(featureCollection([feature(clippedFilter), feature(geo)])); - return inter == null ? 0 : area(inter.geometry); + return polygonsIntersectionArea(geo, filterPolygons); } case "bbox": { const clipped = dropEmptyRings(bboxClip(geo, [spatialFilter.west, spatialFilter.south, spatialFilter.east, spatialFilter.north]).geometry); diff --git a/test/helpers/area.test.ts b/test/helpers/area.test.ts new file mode 100644 index 00000000..55a02c4d --- /dev/null +++ b/test/helpers/area.test.ts @@ -0,0 +1,17 @@ +import { describe, expect, it } from "vitest"; + +import { clipRingToTriangle, sphericalRingArea } from "../../src/helpers/area.js"; + +describe("helpers/area", () => { + describe("clipRingToTriangle", () => { + it("should clip by a triangle whatever its winding", () => { + // The triangle covers the north-eastern quarter of the square. + const square = [[0, 0], [2, 0], [2, 2], [0, 2], [0, 0]]; + const quarter = sphericalRingArea([[1, 1], [2, 1], [2, 2], [1, 2], [1, 1]]); + const [a, b, c] = [[1, 1], [3, 1], [1, 3]]; + + expect(sphericalRingArea(clipRingToTriangle(square, [a, b, c]))).toBeCloseTo(quarter, 12); + expect(sphericalRingArea(clipRingToTriangle(square, [a, c, b]))).toBeCloseTo(quarter, 12); + }); + }); +}); diff --git a/test/wfs/response.test.ts b/test/wfs/response.test.ts index 0c5a01a6..05f4449f 100644 --- a/test/wfs/response.test.ts +++ b/test/wfs/response.test.ts @@ -1,5 +1,9 @@ import { performance } from "node:perf_hooks"; import { describe, expect, it } from "vitest"; +import { intersect } from "@turf/intersect"; +import { feature as turfFeature, featureCollection } from "@turf/helpers"; +import type { Geometry, MultiPolygon, Polygon } from "geojson"; +import area from "../../src/helpers/area"; import { mapToFlatItems, @@ -8,7 +12,6 @@ import { postProcessFeatureCollection, } from "../../src/wfs/response"; import { type SpatialExtraOptions } from "../../src/wfs/schema"; -import type { Geometry } from "geojson"; describe("wfs_engine/response", () => { function getFeatures( @@ -149,7 +152,7 @@ describe("wfs_engine/response", () => { const polygonFeatures = getFeatures(polygonResult); expect(polygonFeatures[0].length).toBeNull(); - expect(polygonFeatures[0].area as number).toBeCloseTo(81361416.69722056, 6); + expect(polygonFeatures[0].area as number).toBeCloseTo(81361415.96674393, 6); }); it("should compute non-zero distance_to_filter_center and intersection_area when a spatial filter is provided", () => { @@ -209,41 +212,47 @@ describe("wfs_engine/response", () => { expect(features[0].distance_to_filter_center as number).toBeCloseTo(2340.9971606708805, 6); }); - it("should stay fast when computing centroid, area, and intersection_area for a large region against many polygons", () => { + it("should stay fast when computing centroid, area, and intersection_area for a large region against many polygons crossing its boundary", () => { const regionCenterLon = 2.35; const regionCenterLat = 48.85; const regionVertices = 3000; const polygonSides = 20; - const featureCount = 100; + const featureCount = 5000; + + /** Close a ring by repeating its first position, as GeoJSON requires. */ + function closeRing(positions: T[]) { + return [...positions, positions[0]]; + } function regularRing(centerLon: number, centerLat: number, radius: number, sides: number) { - return Array.from({ length: sides }, (_, index) => { + return closeRing(Array.from({ length: sides }, (_, index) => { const angle = (Math.PI * 2 * index) / sides - Math.PI / 2; return [centerLon + radius * Math.cos(angle), centerLat + radius * Math.sin(angle)] as [number, number]; - }); + })); + } + + /** Position on the region boundary at the given angle from its center. */ + function regionBoundary(angle: number) { + const radius = 0.22 + 0.03 * Math.sin(7 * angle); + return [ + regionCenterLon + radius * Math.cos(angle), + regionCenterLat + radius * Math.sin(angle), + ]; } const regionGeometry = { type: "Polygon" as const, coordinates: [ - Array.from({ length: regionVertices }, (_, index) => { - const angle = (Math.PI * 2 * index) / regionVertices; - const radius = 0.22 + 0.03 * Math.sin(7 * angle); - return [ - regionCenterLon + radius * Math.cos(angle), - regionCenterLat + radius * Math.sin(angle), - ]; - }), + closeRing(Array.from({ length: regionVertices }, (_, index) => regionBoundary((Math.PI * 2 * index) / regionVertices))), ], }; + // Every polygon is centered on the region boundary, so that each one needs + // an actual intersection to be computed. const featureCollection = { type: "FeatureCollection", features: Array.from({ length: featureCount }, (_, index) => { - const column = index % 100; - const row = Math.floor(index / 100) % 50; - const centerLon = regionCenterLon + (column - 49.5) * 0.0035; - const centerLat = regionCenterLat + (row - 24.5) * 0.0028; + const [centerLon, centerLat] = regionBoundary((Math.PI * 2 * index) / featureCount); const radius = 0.002 + (index % 7) * 0.00025; return { @@ -273,10 +282,13 @@ describe("wfs_engine/response", () => { const features = getFeatures(result); expect(features).toHaveLength(featureCount); expect(features[0].centroid).toMatchObject({ lon: expect.any(Number), lat: expect.any(Number) }); - expect(features[0].area).toEqual(expect.any(Number)); - expect(features[0].intersection_area).toEqual(expect.any(Number)); - expect(features[0].intersection_area).toBeGreaterThanOrEqual(0); - expect(elapsedMs).toBeLessThan(500); + const notPartiallyIntersecting = features.filter((feature) => !( + typeof feature.intersection_area === "number" + && feature.intersection_area > 0 + && feature.intersection_area < (feature.area as number) + )); + expect(notPartiallyIntersecting).toEqual([]); + expect(elapsedMs).toBeLessThan(1000); }); // Regression: `dwithin_point` used to short-circuit `intersection_area` and @@ -320,16 +332,15 @@ describe("wfs_engine/response", () => { }); it("should return the geometry's own area when it is fully inside the filter disc", () => { + // Not a rectangle, whose triangles happen to add up exactly to its area. const containedPolygon = { type: "Polygon", - coordinates: [[ - [1.9999, 47.9999], [2.0001, 47.9999], [2.0001, 48.0001], [1.9999, 48.0001], [1.9999, 47.9999], - ]], + coordinates: [[[2.0003, 48], [2.0001, 48.0003], [1.9997, 48.0002], [1.9997, 47.9998], [2.0001, 47.9997], [2.0003, 48]]], }; const feature = deriveExtras(containedPolygon, ["area", "intersection_area"]); - expect(feature.intersection_area as number).toBeCloseTo(feature.area as number, 6); + expect(feature.intersection_area).toEqual(feature.area); expect(feature.intersection_area as number).toBeLessThanOrEqual(DISC_AREA); }); @@ -471,7 +482,6 @@ describe("wfs_engine/response", () => { expect(feature.intersection_area).toEqual(0); }); - it("should sum the linear parts of a GeometryCollection for length and ignore the other parts", () => { const line = { type: "LineString", coordinates: [[2.3, 48.8], [2.31, 48.81]] }; const lineLength = derive(line, ["length"]).length as number; @@ -511,6 +521,167 @@ describe("wfs_engine/response", () => { expect(feature.length).toBeNull(); }); }); + + describe("intersection_area with an intersects_feature filter", () => { + function rectangle(west: number, south: number, east: number, north: number) { + return [[west, south], [east, south], [east, north], [west, north], [west, south]]; + } + + const referenceSquare = { type: "Polygon" as const, coordinates: [rectangle(2, 48, 3, 49)] }; + + const star: Polygon = { + type: "Polygon", + coordinates: [[ + [2.59, 48.5], [2.5697, 48.5071], [2.5693, 48.5285], [2.5563, 48.5114], [2.5357, 48.5176], + [2.548, 48.5], [2.5357, 48.4824], [2.5563, 48.4886], [2.5693, 48.4715], [2.5697, 48.4929], [2.59, 48.5], + ]], + }; + + function deriveExtras(geometry: unknown, referenceGeometry: Geometry = referenceSquare) { + const result = transformFeatureCollectionResponse({ + type: "FeatureCollection", + features: [{ id: "poly.1", geometry, properties: { name: "poly" } }], + }, { + typename: "TEST:type", + spatial_extras: ["area", "intersection_area"], + intersects_feature_filter: { typename: "TEST:region", feature_id: "region.1" }, + } as Parameters[1], referenceGeometry); + + return getFeatures(result)[0] as { area: number, intersection_area: number }; + } + + it("should return the geometry's own area when it is fully inside the reference", () => { + // Not a rectangle, whose triangles happen to add up exactly to its area. + const feature = deriveExtras(star); + + expect(feature.intersection_area).toEqual(feature.area); + }); + + it("should clip a geometry crossing the reference boundary", () => { + // Split in its middle by the eastern edge of the reference, which is a meridian. + const feature = deriveExtras({ type: "Polygon", coordinates: [rectangle(2.9, 48.4, 3.1, 48.6)] }); + + expect(feature.intersection_area / feature.area).toBeCloseTo(0.5, 9); + }); + + it("should return the reference's area for a geometry containing it", () => { + const feature = deriveExtras({ type: "Polygon", coordinates: [rectangle(1, 47, 4, 50)] }); + + expect(feature.intersection_area).toBeCloseTo(area(referenceSquare), 3); + }); + + it("should exclude the holes of the reference", () => { + const hole = rectangle(2.4, 48.4, 2.6, 48.6); + const referenceWithHole = { type: "Polygon" as const, coordinates: [rectangle(2, 48, 3, 49), hole] }; + const feature = deriveExtras({ type: "Polygon", coordinates: [rectangle(2.3, 48.3, 2.7, 48.7)] }, referenceWithHole); + + const holeArea = area({ type: "Polygon", coordinates: [hole] }); + expect(feature.intersection_area).toBeCloseTo(feature.area - holeArea, 3); + }); + + it("should sum the intersections with every part of a MultiPolygon reference", () => { + const referenceMultiPolygon = { + type: "MultiPolygon" as const, + coordinates: [[rectangle(2, 48, 3, 49)], [rectangle(3.2, 48, 4, 49)]], + }; + // Overlaps [2.9, 3] with the first part and [3.2, 3.3] with the second. + const feature = deriveExtras({ type: "Polygon", coordinates: [rectangle(2.9, 48.4, 3.3, 48.6)] }, referenceMultiPolygon); + + expect(feature.intersection_area / feature.area).toBeCloseTo(0.5, 9); + }); + + it("should sum the intersections of every part of a MultiPolygon geometry", () => { + // One part only touching the eastern edge of the reference, the other inside it. + const inside = rectangle(2.4, 48.4, 2.6, 48.6); + const feature = deriveExtras({ type: "MultiPolygon", coordinates: [[rectangle(3, 48.4, 3.2, 48.6)], [inside]] }); + + expect(feature.intersection_area).toEqual(area({ type: "Polygon", coordinates: [inside] })); + }); + + /** Area of the intersection computed by general polygon clipping, as a reference. */ + function polyclipIntersectionArea(geometry: Polygon | MultiPolygon, referenceGeometry: Polygon | MultiPolygon) { + const inter = intersect(featureCollection([turfFeature(geometry), turfFeature(referenceGeometry)])); + return inter == null ? 0 : area(inter.geometry); + } + + it("should clip a non-convex geometry crossing the reference boundary", () => { + // U shape whose both branches cross the eastern edge of the reference. + const uShape: Polygon = { + type: "Polygon", + coordinates: [[[2.8, 48.4], [3.2, 48.4], [3.2, 48.45], [2.9, 48.45], [2.9, 48.55], [3.2, 48.55], [3.2, 48.6], [2.8, 48.6], [2.8, 48.4]]], + }; + const feature = deriveExtras(uShape); + + expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(uShape, referenceSquare), 0); + }); + + it("should exclude the holes of the geometry", () => { + const withHole: Polygon = { + type: "Polygon", + coordinates: [rectangle(2.8, 48.4, 3.2, 48.6), rectangle(2.9, 48.45, 3.1, 48.55)], + }; + const feature = deriveExtras(withHole); + + expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(withHole, referenceSquare), 0); + }); + + it("should fall back on general polygon clipping for a self-intersecting geometry", () => { + // Bow tie, which cannot be triangulated. + const bowTie: Polygon = { + type: "Polygon", + coordinates: [[[2.8, 48.4], [3.2, 48.6], [3.2, 48.4], [2.8, 48.6], [2.8, 48.4]]], + }; + const feature = deriveExtras(bowTie); + + expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(bowTie, referenceSquare), 0); + + // A reference in the empty wedge below the crossing point is not intersected. + const inWedge = { type: "Polygon" as const, coordinates: [rectangle(2.98, 48.41, 3.02, 48.43)] }; + expect(deriveExtras(bowTie, inWedge).intersection_area).toEqual(0); + }); + + it("should return null for a non-areal geometry", () => { + const feature = deriveExtras({ type: "LineString", coordinates: [[2.5, 48.5], [3.5, 48.5]] }); + + expect(feature.intersection_area).toBeNull(); + }); + + it("should clip a non-convex geometry crossing a non-convex reference boundary", () => { + // Notch on the eastern side of the reference, whose tip lies inside the star. + const notchedReference: Polygon = { + type: "Polygon", + coordinates: [[[2.4, 48.4], [2.6, 48.4], [2.55, 48.5], [2.6, 48.6], [2.4, 48.6], [2.4, 48.4]]], + }; + const feature = deriveExtras(star, notchedReference); + + expect(feature.intersection_area / feature.area).toBeCloseTo(0.2856, 4); + expect(Math.abs(feature.intersection_area - polyclipIntersectionArea(star, notchedReference))).toBeLessThanOrEqual(1e-6 * feature.area); + }); + + it("should stay precise on long thin geometries crossing a winding reference boundary", () => { + // Clipping leaves zero-width spikes along the long triangles of such geometries, + // which must not contribute to the area. + const winding = Array.from({ length: 500 }, (_, index) => { + const angle = (Math.PI * 2 * index) / 500; + const radius = 0.2 * (1 + 0.05 * Math.sin(53 * angle)); + return [2.5 + radius * Math.cos(angle), 48.5 + radius * Math.sin(angle)]; + }); + const reference: Polygon = { type: "Polygon", coordinates: [[...winding, winding[0]]] }; + + // 10 km long and 20 m wide diagonal strip, with a vertex every 100 m, centered on the boundary. + const side = (offset: number) => Array.from({ length: 101 }, (_, index) => { + const along = -0.032 + (0.064 * index) / 100; + return [2.7 + along - offset, 48.5 + along + offset]; + }); + const ring = [...side(0.00007), ...side(-0.00007).reverse()]; + const strip: Polygon = { type: "Polygon", coordinates: [[...ring, ring[0]]] }; + + const feature = deriveExtras(strip, reference); + + const expected = polyclipIntersectionArea(strip, reference); + expect(Math.abs(feature.intersection_area - expected)).toBeLessThanOrEqual(1e-6 * feature.area); + }); + }); }); // --- postProcessFeatureCollection --- From 51483171320796031d0f07f65fa9fecfb1191f18 Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Tue, 29 Sep 2026 16:50:23 +0200 Subject: [PATCH 2/6] perf: add FilterClipper for intersection_area computation --- src/wfs/spatialExtras.ts | 63 +++++++++++++++++++++++++++++++-------- test/wfs/response.test.ts | 9 ++++++ 2 files changed, 59 insertions(+), 13 deletions(-) diff --git a/src/wfs/spatialExtras.ts b/src/wfs/spatialExtras.ts index a490fc1a..fb190128 100644 --- a/src/wfs/spatialExtras.ts +++ b/src/wfs/spatialExtras.ts @@ -4,7 +4,7 @@ import turfLength from "@turf/length"; import { intersect } from "@turf/intersect"; import { circle } from "@turf/circle"; import { bboxPolygon } from "@turf/bbox-polygon"; -import type { Geometry, LineString, MultiLineString, MultiPolygon, Point, Polygon, Position } from "geojson"; +import type { BBox, Geometry, LineString, MultiLineString, MultiPolygon, Point, Polygon, Position } from "geojson"; import distance from "../helpers/distance.js"; import { feature, featureCollection } from "@turf/helpers"; import { getSpatialFilter } from "./spatialFilter.js"; @@ -134,7 +134,7 @@ function isComputableGeometry(geometry: unknown) : geometry is Geometry { type SpatialContext = { filterCentroid: Point | null, - filterPolygons: Polygon | MultiPolygon | null + clipFilter: FilterClipper | null } export function prepareSpatialContext(input: FeatureCollectionPostProcessInput, resolvedGeometryRef?: Geometry) : SpatialContext { @@ -143,7 +143,7 @@ export function prepareSpatialContext(input: FeatureCollectionPostProcessInput, const context : SpatialContext = { filterCentroid: null, - filterPolygons: null, + clipFilter: null, }; if (!requires_distance_to_filter_center && !requires_intersection_area) { @@ -164,7 +164,8 @@ export function prepareSpatialContext(input: FeatureCollectionPostProcessInput, if (requires_intersection_area) { try { - context.filterPolygons = geometryToPolygons(spatialFilterToGeometry(spatialFilter, resolvedGeometryRef)); + const filterPolygons = geometryToPolygons(spatialFilterToGeometry(spatialFilter, resolvedGeometryRef)); + context.clipFilter = filterPolygons && makeFilterClipper(filterPolygons); } catch {} } @@ -190,19 +191,52 @@ export function dropEmptyRings(geom: Geometry) : Polygon | MultiPolygon | null { return null; } +/** Return the filter clipped to a tile containing the bbox, or null if nothing areal remains. */ +type FilterClipper = (bbox: BBox) => Polygon | MultiPolygon | null; + +/** Build a `FilterClipper` which only processes the filter vertices near the bbox. + * + * The filter is clipped to nested tiles, each from the cached clip of its parent. + * Tiles are 1.5 times as wide as their cell, so that a bbox at most half a cell + * wide lies in the tile of the cell containing its south-west corner. + */ +function makeFilterClipper(filter: Polygon | MultiPolygon) : FilterClipper { + const [west, south, east, north] = bbox(filter); + const size = Math.max(east - west, north - south); + const tiles = new Map(); + function tile(zoom: number, x: number, y: number) : Polygon | MultiPolygon | null { + if (zoom == 0) return filter; + const key = `${zoom}/${x}/${y}`; + if (!tiles.has(key)) { + const parent = tile(zoom - 1, Math.floor(x / 2), Math.floor(y / 2)); + const cell = size / 2 ** zoom; + tiles.set(key, parent && dropEmptyRings(bboxClip(parent, [west + x * cell, south + y * cell, west + (x + 1.5) * cell, south + (y + 1.5) * cell]).geometry)); + } + return tiles.get(key) ?? null; + } + return ([minX, minY, maxX, maxY]) => { + // Deepest tile whose cell is at least twice as wide as the bbox. + const width = Math.max(maxX - minX, maxY - minY); + let zoom = 0; + while (zoom < 20 && 2 * width <= size / 2 ** (zoom + 1)) zoom++; + const cell = size / 2 ** zoom; + return tile(zoom, Math.floor((minX - west) / cell), Math.floor((minY - south) / cell)); + }; +} /** Return the area of the intersection between two areal (2D) geometries. * * Clipping by a triangle is simple, fast and robust, unlike general polygon - * clipping: `geo` is split into triangles, by which `filterPolygons` is clipped. + * clipping: `geo` is split into triangles, by which the filter is clipped. */ -function polygonsIntersectionArea(geo: Polygon | MultiPolygon, filterPolygons: Polygon | MultiPolygon) : number { +function polygonsIntersectionArea(geo: Polygon | MultiPolygon, clipFilter: FilterClipper) : number { // Whatever lies outside the feature's bbox cannot intersect it, so clipping the // reference first leaves the result unchanged while only the neighbouring // vertices are processed afterwards. - const clippedFilter = dropEmptyRings(bboxClip(filterPolygons, bbox(geo)).geometry); + const geoBbox = bbox(geo); + const nearFilter = clipFilter(geoBbox); + const clippedFilter = nearFilter && dropEmptyRings(bboxClip(nearFilter, geoBbox).geometry); if (!clippedFilter) return 0; // no areal overlap - const filterParts = clippedFilter.type == "Polygon" ? [clippedFilter.coordinates] : clippedFilter.coordinates; let total = 0; for (const polygon of geo.type == "Polygon" ? [geo.coordinates] : geo.coordinates) { @@ -216,7 +250,10 @@ function polygonsIntersectionArea(geo: Polygon | MultiPolygon, filterPolygons: P } let covered = 0; for (const triangle of triangles) { - for (const [outer, ...holes] of filterParts) { + // For a large polygon, clipping the filter near each triangle is faster. + const nearTriangle = triangles.length >= 128 ? clippedFilter : clipFilter(bbox({ type: "MultiPoint", coordinates: triangle })); + if (!nearTriangle) continue; + for (const [outer, ...holes] of nearTriangle.type == "Polygon" ? [nearTriangle.coordinates] : nearTriangle.coordinates) { covered += area({ type: "Polygon", coordinates: [clipRingToTriangle(outer, triangle)]}); for (const hole of holes) { covered -= area({ type: "Polygon", coordinates: [clipRingToTriangle(hole, triangle)]}); @@ -236,7 +273,7 @@ function polygonsIntersectionArea(geo: Polygon | MultiPolygon, filterPolygons: P * or the filter geometry could not be prepared. 0 only when both are areal and * do not overlap. */ -function intersectionAreaWithSpatialFilter(geom: Geometry, spatialFilter: SpatialFilter, filterPolygons: Polygon | MultiPolygon | null) : number | null { +function intersectionAreaWithSpatialFilter(geom: Geometry, spatialFilter: SpatialFilter, clipFilter: FilterClipper | null) : number | null { const geo = geometryToPolygons(geom); if (!geo) { return null; // non-areal feature @@ -249,8 +286,8 @@ function intersectionAreaWithSpatialFilter(geom: Geometry, spatialFilter: Spatia case "dwithin_point": case "intersects_feature": case "travel_time": { - if (!filterPolygons) return null; // non-areal filter, or filter preparation failed - return polygonsIntersectionArea(geo, filterPolygons); + if (!clipFilter) return null; // non-areal filter, or filter preparation failed + return polygonsIntersectionArea(geo, clipFilter); } case "bbox": { const clipped = dropEmptyRings(bboxClip(geo, [spatialFilter.west, spatialFilter.south, spatialFilter.east, spatialFilter.north]).geometry); @@ -353,7 +390,7 @@ export function deriveFromGeometry(geometry: unknown, input: FeatureCollectionPo if (requires_intersection_area) { try { - ret.intersection_area = intersectionAreaWithSpatialFilter(geo, spatialFilter, context.filterPolygons); + ret.intersection_area = intersectionAreaWithSpatialFilter(geo, spatialFilter, context.clipFilter); } catch { ret.intersection_area = null; } diff --git a/test/wfs/response.test.ts b/test/wfs/response.test.ts index 05f4449f..b3fa169b 100644 --- a/test/wfs/response.test.ts +++ b/test/wfs/response.test.ts @@ -681,6 +681,15 @@ describe("wfs_engine/response", () => { const expected = polyclipIntersectionArea(strip, reference); expect(Math.abs(feature.intersection_area - expected)).toBeLessThanOrEqual(1e-6 * feature.area); }); + + it("should handle geometries smaller than the deepest tile of the reference", () => { + // About 1 cm wide, inside the reference and in an empty tile outside it. + const tinySquare = (lon: number, lat: number) => ({ type: "Polygon", coordinates: [rectangle(lon, lat, lon + 1e-7, lat + 1e-7)] }); + const inside = deriveExtras(tinySquare(2.5, 48.5)); + + expect(Math.abs(inside.intersection_area - inside.area)).toBeLessThanOrEqual(1e-6 * inside.area); + expect(deriveExtras(tinySquare(3.5, 48.5)).intersection_area).toEqual(0); + }); }); }); From f720d5eca46616c97a03518fbf10bdced0f0de6a Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Wed, 30 Sep 2026 12:23:27 +0200 Subject: [PATCH 3/6] perf: subdivide boxes for intersection_area --- src/wfs/spatialExtras.ts | 104 +++++++++++++++++++++++++++----------- test/wfs/response.test.ts | 59 +++++++++++++++------ 2 files changed, 119 insertions(+), 44 deletions(-) diff --git a/src/wfs/spatialExtras.ts b/src/wfs/spatialExtras.ts index fb190128..0fec7650 100644 --- a/src/wfs/spatialExtras.ts +++ b/src/wfs/spatialExtras.ts @@ -224,44 +224,90 @@ function makeFilterClipper(filter: Polygon | MultiPolygon) : FilterClipper { }; } -/** Return the area of the intersection between two areal (2D) geometries. +/** Polygons of a polygonal geometry, as lists of rings. */ +function polygonsOf(geo: Polygon | MultiPolygon) : Position[][][] { + return geo.type == "Polygon" ? [geo.coordinates] : geo.coordinates; +} + +function positionCount(geo: Polygon | MultiPolygon) : number { + return polygonsOf(geo).flat().reduce((count, ring) => count + ring.length, 0); +} + +/** True when `geo` is the whole `box`, as clipping returns a geometry covering it. */ +function isBox(geo: Polygon | MultiPolygon, [west, south, east, north]: BBox) : boolean { + const rings = polygonsOf(geo).flat(); + if (rings.length != 1 || rings[0].length != 5) return false; + const corners = rings[0].slice(0, -1); + return corners.every(([x, y]) => (x == west || x == east) && (y == south || y == north)) + && new Set(corners.map(String)).size == 4; +} + +/** Beyond this many positions, geometries are split before being triangulated. */ +const MAX_TRIANGULATED_POSITIONS = 32; +/** Bound on the number of halvings, in case positions pile up at the same place. */ +const MAX_HALVINGS = 24; + +/** Return the area of the intersection between two areal geometries lying in `box`. * * Clipping by a triangle is simple, fast and robust, unlike general polygon - * clipping: `geo` is split into triangles, by which the filter is clipped. + * clipping: the geometry with fewer positions is split into triangles, by which + * the other one is clipped. As this costs the product of their numbers of + * positions, the box is first halved until one of them is small. */ -function polygonsIntersectionArea(geo: Polygon | MultiPolygon, clipFilter: FilterClipper) : number { - // Whatever lies outside the feature's bbox cannot intersect it, so clipping the - // reference first leaves the result unchanged while only the neighbouring - // vertices are processed afterwards. - const geoBbox = bbox(geo); - const nearFilter = clipFilter(geoBbox); - const clippedFilter = nearFilter && dropEmptyRings(bboxClip(nearFilter, geoBbox).geometry); - if (!clippedFilter) return 0; // no areal overlap +function boxIntersectionArea(a: Polygon | MultiPolygon, b: Polygon | MultiPolygon, box: BBox, halvings = 0) : number { + const [small, large] = positionCount(a) <= positionCount(b) ? [a, b] : [b, a]; + if (positionCount(small) > MAX_TRIANGULATED_POSITIONS && halvings < MAX_HALVINGS) { + const [west, south, east, north] = box; + const halves : BBox[] = east - west > north - south + ? [[west, south, (west + east) / 2, north], [(west + east) / 2, south, east, north]] + : [[west, south, east, (south + north) / 2], [west, (south + north) / 2, east, north]]; + return halves.reduce((total, half) => { + const smallHalf = dropEmptyRings(bboxClip(small, half).geometry); + const largeHalf = smallHalf && dropEmptyRings(bboxClip(large, half).geometry); + return largeHalf ? total + boxIntersectionArea(smallHalf, largeHalf, half, halvings + 1) : total; + }, 0); + } + + if (isBox(small, box)) return area(large); + const triangulations = polygonsOf(small).map(triangulate); + if (triangulations.includes(null)) { + // Fall back on general polygon clipping, as for self-intersecting rings. + const inter = intersect(featureCollection([feature(small), feature(large)])); + return inter == null ? 0 : area(inter.geometry); + } + let total = 0; + for (const triangle of triangulations.flatMap((triangles) => triangles!)) { + total += area({ type: "MultiPolygon", coordinates: polygonsOf(large).map((rings) => rings.map((ring) => clipRingToTriangle(ring, triangle))) }); + } + return total; +} +/** Return the area of the intersection between two areal (2D) geometries. */ +function polygonsIntersectionArea(geo: Polygon | MultiPolygon, clipFilter: FilterClipper) : number { let total = 0; - for (const polygon of geo.type == "Polygon" ? [geo.coordinates] : geo.coordinates) { + for (const polygon of polygonsOf(geo)) { const polygonGeometry : Polygon = { type: "Polygon", coordinates: polygon }; - const triangles = triangulate(polygon); - if (!triangles) { - // Fall back on general polygon clipping. - const inter = intersect(featureCollection([feature(clippedFilter), feature(polygonGeometry)])); - total += inter == null ? 0 : area(inter.geometry); + // Whatever lies outside the bbox of one geometry cannot intersect it, so + // clipping the other first leaves the result unchanged while only the + // neighbouring positions are processed afterwards. + const polygonBbox = bbox(polygonGeometry); + const nearFilter = clipFilter(polygonBbox); + const clippedFilter = nearFilter && dropEmptyRings(bboxClip(nearFilter, polygonBbox).geometry); + if (!clippedFilter) continue; // no areal overlap + const polygonArea = area(polygonGeometry); + if (isBox(clippedFilter, polygonBbox)) { // the filter contains the polygon + total += polygonArea; continue; } - let covered = 0; - for (const triangle of triangles) { - // For a large polygon, clipping the filter near each triangle is faster. - const nearTriangle = triangles.length >= 128 ? clippedFilter : clipFilter(bbox({ type: "MultiPoint", coordinates: triangle })); - if (!nearTriangle) continue; - for (const [outer, ...holes] of nearTriangle.type == "Polygon" ? [nearTriangle.coordinates] : nearTriangle.coordinates) { - covered += area({ type: "Polygon", coordinates: [clipRingToTriangle(outer, triangle)]}); - for (const hole of holes) { - covered -= area({ type: "Polygon", coordinates: [clipRingToTriangle(hole, triangle)]}); - } - } - } + const filterBbox = bbox(clippedFilter); + // Clipping the polygon to its own bbox would only copy it. + const clippedPolygon = filterBbox.every((value, index) => value == polygonBbox[index]) + ? polygonGeometry + : dropEmptyRings(bboxClip(polygonGeometry, filterBbox).geometry); + if (!clippedPolygon) continue; + + const covered = boxIntersectionArea(clippedPolygon, clippedFilter, filterBbox); // Snap to 0 or to the whole polygon despite rounding errors. - const polygonArea = area(polygonGeometry); total += covered > polygonArea * (1 - 1e-9) ? polygonArea : covered > polygonArea * 1e-9 ? covered : 0; } return total; diff --git a/test/wfs/response.test.ts b/test/wfs/response.test.ts index b3fa169b..cb2537c9 100644 --- a/test/wfs/response.test.ts +++ b/test/wfs/response.test.ts @@ -291,6 +291,42 @@ describe("wfs_engine/response", () => { expect(elapsedMs).toBeLessThan(1000); }); + it("should stay fast when computing intersection_area for a large geometry with holes crossing a detailed reference boundary", () => { + // Like a forest with clearings across the boundary of an isochrone: both + // geometries have many positions in the same place, which used to cost + // their product (over a second here). + function wavyRing(centerLon: number, centerLat: number, radius: number, vertices: number, waves: number) { + const positions = Array.from({ length: vertices }, (_, index) => { + const angle = (Math.PI * 2 * index) / vertices; + const wavyRadius = radius * (1 + 0.02 * Math.sin(waves * angle)); + return [centerLon + wavyRadius * Math.cos(angle), centerLat + wavyRadius * Math.sin(angle)]; + }); + return [...positions, positions[0]]; + } + + const reference: Polygon = { type: "Polygon", coordinates: [wavyRing(2.35, 48.85, 0.3, 20000, 200)] }; + // Centered on the eastern boundary of the reference, with clearings on a grid. + const clearings = Array.from({ length: 100 }, (_, index) => [2.65 + 0.014 * (index % 10 - 4.5), 48.85 + 0.014 * (Math.floor(index / 10) - 4.5)]) + .filter(([lon, lat]) => Math.hypot(lon - 2.65, lat - 48.85) < 0.08) + .map(([lon, lat]) => wavyRing(lon, lat, 0.003, 16, 3).reverse()); + const forest: Polygon = { type: "Polygon", coordinates: [wavyRing(2.65, 48.85, 0.1, 20000, 300), ...clearings] }; + + const start = performance.now(); + const result = transformFeatureCollectionResponse({ + type: "FeatureCollection", + features: [{ id: "forest.1", geometry: forest, properties: {} }], + }, { + typename: "TEST:type", + spatial_extras: ["intersection_area"], + intersects_feature_filter: { typename: "TEST:region", feature_id: "region.1" }, + }, reference); + const elapsedMs = performance.now() - start; + + const inter = intersect(featureCollection([turfFeature(forest), turfFeature(reference)])); + expect(getFeatures(result)[0].intersection_area as number).toBeCloseTo(area(inter!.geometry), 3); + expect(elapsedMs).toBeLessThan(300); + }); + // Regression: `dwithin_point` used to short-circuit `intersection_area` and // return the feature's own geometry, on the false premise that a DWITHIN // match implies containment. DWITHIN matches as soon as ANY part of the @@ -557,6 +593,14 @@ describe("wfs_engine/response", () => { expect(feature.intersection_area).toEqual(feature.area); }); + it("should return the geometry's own area when it is inside the reference but not its bbox", () => { + // The hypotenuse of the reference cuts the north-eastern corner of the bbox of the star, not the star. + const reference = { type: "Polygon" as const, coordinates: [[[2, 48], [2.97, 48], [2, 49.3115], [2, 48]]] }; + const feature = deriveExtras(star, reference); + + expect(feature.intersection_area).toEqual(feature.area); + }); + it("should clip a geometry crossing the reference boundary", () => { // Split in its middle by the eastern edge of the reference, which is a meridian. const feature = deriveExtras({ type: "Polygon", coordinates: [rectangle(2.9, 48.4, 3.1, 48.6)] }); @@ -625,21 +669,6 @@ describe("wfs_engine/response", () => { expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(withHole, referenceSquare), 0); }); - it("should fall back on general polygon clipping for a self-intersecting geometry", () => { - // Bow tie, which cannot be triangulated. - const bowTie: Polygon = { - type: "Polygon", - coordinates: [[[2.8, 48.4], [3.2, 48.6], [3.2, 48.4], [2.8, 48.6], [2.8, 48.4]]], - }; - const feature = deriveExtras(bowTie); - - expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(bowTie, referenceSquare), 0); - - // A reference in the empty wedge below the crossing point is not intersected. - const inWedge = { type: "Polygon" as const, coordinates: [rectangle(2.98, 48.41, 3.02, 48.43)] }; - expect(deriveExtras(bowTie, inWedge).intersection_area).toEqual(0); - }); - it("should return null for a non-areal geometry", () => { const feature = deriveExtras({ type: "LineString", coordinates: [[2.5, 48.5], [3.5, 48.5]] }); From fe7a2f2a7ce19c431b9e85d3c10e6aebb9455554 Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:38:22 +0200 Subject: [PATCH 4/6] perf(intersection_area): integrate both geometries over a recursive box split Replace the per-feature triangle clipping (FilterClipper and box halving) with src/helpers/intersectionArea.ts: - the filter rings are copied once into flat arrays; the tile cache clips them with an exact axis-aligned Sutherland-Hodgman instead of bboxClip - the feature is clipped to the bbox of the filter near it, and the common box is split in halves until one side has few positions; a box that one side covers entirely, or not at all, ends the recursion - at the leaves, the smaller side is triangulated and the other side's rings are clipped by each triangle; areas come from winding numbers, exact on the sphere and computed relative to a local origin - self-intersecting rings and oversized leaves fall back on @turf/intersect, which restores the bow tie test removed by 3615e1c Median time per request on real Geoplateforme data, 3615e1c -> this commit: 2.1 s -> 0.47 s (communes x departement 35), 1.9 s -> 0.36 s (departements x region 24), 636 -> 223 ms (5000 vegetation zones x departement 35). Results match GEOS 3.13 within 1e-8 of the feature area. - move dropEmptyRings to helpers/geojson.ts - drop clipRingToTriangle and triangulate from helpers/area.ts, now unused --- src/helpers/area.ts | 58 +- src/helpers/geojson.ts | 21 +- src/helpers/intersectionArea.ts | 793 ++++++++++++++++++++++++++ src/wfs/spatialExtras.ts | 164 +----- test/helpers/area.test.ts | 23 +- test/helpers/intersectionArea.test.ts | 69 +++ test/wfs/geometry.test.ts | 2 +- test/wfs/response.test.ts | 22 +- 8 files changed, 931 insertions(+), 221 deletions(-) create mode 100644 src/helpers/intersectionArea.ts create mode 100644 test/helpers/intersectionArea.test.ts diff --git a/src/helpers/area.ts b/src/helpers/area.ts index 56ee1dad..c7cb435e 100644 --- a/src/helpers/area.ts +++ b/src/helpers/area.ts @@ -1,63 +1,7 @@ -import earcut, { deviation, flatten } from "earcut"; import type { Geometry, Position } from "geojson"; -type Triangle = [Position, Position, Position]; - const DEGREES_TO_RADIANS = Math.PI / 180; -const EARTH_RADIUS = 6371008.7714; - -/** Positive when `c` lies on the left of the line going from `a` to `b`, negative on its right. */ -function cross(a: Position, b: Position, c: Position) : number { - return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]); -} - -/** Clip a closed ring by a triangle, using the Sutherland–Hodgman algorithm. - * - * The result may contain zero-width spikes along the triangle edges, which do not - * change its area. It is empty when nothing areal survives. - */ -export function clipRingToTriangle(ring: Position[], [a, b, c]: Triangle) : Position[] { - // Orient the triangle counterclockwise, so that its interior lies on the left of each edge. - const edges = cross(a, b, c) > 0 ? [[a, b], [b, c], [c, a]] : [[a, c], [c, b], [b, a]]; - let clipped = ring.slice(0, -1); // Drop the closing position - for (const [start, end] of edges) { - const input = clipped; - clipped = []; - input.forEach((current, index) => { - const previous = input.at(index - 1)!; // The last position for the first one - const previousSide = cross(start, end, previous); - const currentSide = cross(start, end, current); - if ((previousSide >= 0) !== (currentSide >= 0)) { - // The edge crosses the line: keep the crossing point. Both sides have - // different signs, so the denominator cannot be 0. - const t = previousSide / (previousSide - currentSide); - clipped.push([previous[0] + t * (current[0] - previous[0]), previous[1] + t * (current[1] - previous[1])]); - } - if (currentSide >= 0) { - clipped.push(current); - } - }); - } - return clipped.length < 3 ? [] : [...clipped, clipped[0]]; -} - -/** Split a polygon into triangles, or return null if it cannot be done reliably, - * as for self-intersecting rings. - */ -export function triangulate(polygon: Position[][]) : Triangle[] | null { - const { vertices, holes, dimensions } = flatten(polygon); - const indices = earcut(vertices, holes, dimensions); - // The triangles must cover the polygon exactly. - if (deviation(vertices, holes, dimensions, indices) > 1e-6) { - return null; - } - const position = (index: number) : Position => [vertices[index * dimensions], vertices[index * dimensions + 1]]; - const triangles : Triangle[] = []; - for (let i = 0; i < indices.length; i += 3) { - triangles.push([position(indices[i]), position(indices[i + 1]), position(indices[i + 2])]); - } - return triangles; -} +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). * diff --git a/src/helpers/geojson.ts b/src/helpers/geojson.ts index da67f1a0..0a7eb179 100644 --- a/src/helpers/geojson.ts +++ b/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 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; +} diff --git a/src/helpers/intersectionArea.ts b/src/helpers/intersectionArea.ts new file mode 100644 index 00000000..2fb270cf --- /dev/null +++ b/src/helpers/intersectionArea.ts @@ -0,0 +1,793 @@ +import earcut, { deviation } from "earcut"; +import { bbox } from "@turf/bbox"; +import { bboxClip } from "@turf/bbox-clip"; +import { intersect } from "@turf/intersect"; +import { feature, featureCollection } from "@turf/helpers"; +import type { BBox, MultiPolygon, Polygon, Position } from "geojson"; +import area, { EARTH_RADIUS } from "./area.js"; +import { dropEmptyRings } from "./geojson.js"; + +/** + * Area of the part of areal features lying inside an areal filter. + * + * The filter is prepared once: its rings are copied into flat arrays and lazily clipped + * to nested tiles, so that a feature only meets the filter positions near its bbox. Both + * geometries are then clipped to their common box, which is split recursively in halves + * until one side has few positions. A box that one side covers entirely, or not at all, + * ends the recursion. At the leaves, the smaller side is triangulated (earcut) and the + * rings of the other side are clipped by each triangle (Sutherland-Hodgman). + * + * Correctness rests on winding numbers: clipping a closed ring by a convex region keeps + * its winding number inside the region and sets it to 0 outside, whatever zero-width + * bridges the clipping leaves, so signed areas stay exact. Each ring carries mu = ±1 (+1 + * for an outer ring, -1 for a hole, times its orientation), the multiplicity of a geometry + * is the sum of mu · winding over its rings, and the intersection area is the integral of + * the product of both multiplicities: the polygons of a MultiPolygon are summed, overlaps + * included, as `area` does. + * + * Self-intersecting rings (a bow tie, a clipped piece with the wrong orientation, a proper + * crossing where earcut fails) and oversized leaves are left to `@turf/intersect`. + * + * Areas are exact on the sphere for edges straight in lon/lat (same integral and radius as + * `area`), evaluated relative to a local origin to avoid cancellation on small rings. + */ + +/** Area (m²) of the part of an areal feature lying inside the prepared filter. */ +export type IntersectionArea = (geo: Polygon | MultiPolygon) => number; + +/** A box becomes a leaf when small × big <= max(LEAF_PRODUCT, LEAF_FACTOR × (small + big)) positions. */ +const LEAF_PRODUCT = 1024; +const LEAF_FACTOR = 16; +/** Bound on the recursion, in case positions pile up at the same place. */ +const MAX_DEPTH = 60; +/** Triangle × position operations of a leaf above which `@turf/intersect` takes over. */ +const MAX_LEAF_OPERATIONS = 1e6; +/** Relative snap to the feature area, and to 0 (tighter, so that real slivers survive). */ +const SNAP_FULL = 1e-9; +const SNAP_ZERO = 1e-12; +/** earcut deviation above which a signed fan is used instead (local coordinates are precise). */ +const MAX_DEVIATION = 1e-9; + +const DEGREES_TO_RADIANS = Math.PI / 180; +const HALF_DEGREES_TO_RADIANS = DEGREES_TO_RADIANS / 2; +const R2 = EARTH_RADIUS ** 2; + +/** Thrown to hand a feature over to `@turf/intersect`. */ +const FALLBACK = Symbol("fallback"); + +/** + * A ring in flat local coordinates `[x0, y0, x1, y1, ...]`, without its closing position, + * with its bbox. `o` is the orientation of the original ring (+1 counterclockwise), `mu` is + * `o` for an outer ring and `-o` for a hole. + */ +type Ring = { c: Float64Array; n: number; x0: number; y0: number; x1: number; y1: number; mu: number; o: number }; + +/** Latitude of the local origin, in radians, with its sine and cosine. */ +type AreaContext = { phi0: number; s0: number; c0: number }; + +type Axis = 0 | 1; + +// Scratch buffers, reused from one clip to the next: no state survives a call. +let bufferA = new Float64Array(1 << 12); +let bufferB = new Float64Array(1 << 12); + +function growA(size: number) : Float64Array { + if (bufferA.length < size) bufferA = new Float64Array(Math.max(size, 2 * bufferA.length)); + return bufferA; +} + +function growB(size: number) : Float64Array { + if (bufferB.length < size) bufferB = new Float64Array(Math.max(size, 2 * bufferB.length)); + return bufferB; +} + +function polygonsOf(geo: Polygon | MultiPolygon) : Position[][][] { + return geo.type == "Polygon" ? [geo.coordinates] : geo.coordinates; +} + +// --- Rings --- + +/** A clipped piece of a ring of orientation `o`, copied from the first `count2` values of `src`. */ +function makeRing(src: Float64Array, count2: number, mu: number, o: number) : Ring { + const c = src.slice(0, count2); + const n = count2 >> 1; + // bbox, and twice the planar signed area as a fan from the first position + const ox = c[0], oy = c[1]; + let x0 = ox, y0 = oy, x1 = ox, y1 = oy; + let area2 = 0, px = 0, py = 0; + for (let i = 2; i < count2; i += 2) { + const x = c[i], y = c[i + 1]; + if (x < x0) x0 = x; + else if (x > x1) x1 = x; + if (y < y0) y0 = y; + else if (y > y1) y1 = y; + const qx = x - ox, qy = y - oy; + area2 += px * qy - qx * py; + px = qx; + py = qy; + } + // The pieces of a simple ring keep its orientation: an opposite one, well beyond + // rounding, reveals a self-intersecting ring. + const extent = Math.max(x1 - x0, y1 - y0); + const tolerance = 2e-6 * (x1 - x0) * (y1 - y0) + 1e-12 * n * extent * extent; + if (area2 * o < -tolerance) throw FALLBACK; + return { c, n, x0, y0, x1, y1, mu, o }; +} + +/** Local copy of a GeoJSON ring (x = lon - lon0, y = lat - lat0) without its closing position, or null below 3 positions. */ +function localRing(positions: Position[], lon0: number, lat0: number) : Ring | null { + let n = positions.length; + if (n > 1) { + const first = positions[0], last = positions[n - 1]; + if (first[0] === last[0] && first[1] === last[1]) n--; + } + if (n < 3) return null; + const c = new Float64Array(2 * n); + let x0 = Infinity, y0 = Infinity, x1 = -Infinity, y1 = -Infinity; + for (let i = 0; i < n; i++) { + const x = positions[i][0] - lon0, y = positions[i][1] - lat0; + c[2 * i] = x; + c[2 * i + 1] = y; + if (x < x0) x0 = x; + if (x > x1) x1 = x; + if (y < y0) y0 = y; + if (y > y1) y1 = y; + } + return { c, n, x0, y0, x1, y1, mu: 0, o: 0 }; +} + +/** Twice the planar signed area (positive counterclockwise), as a fan from the first position. */ +function planarArea2(c: Float64Array, n: number) : number { + const ox = c[0], oy = c[1]; + let s = 0; + let px = c[2] - ox, py = c[3] - oy; + for (let i = 4; i < 2 * n; i += 2) { + const x = c[i] - ox, y = c[i + 1] - oy; + s += px * y - x * py; + px = x; + py = y; + } + return s; +} + +function countPositions(rings: Ring[]) : number { + let n = 0; + for (const r of rings) n += r.n; + return n; +} + +// --- Spherical area in local coordinates --- + +/** + * Signed area on the unit sphere (positive counterclockwise) of a ring in local coordinates, + * for edges straight in lon/lat: the sum of dLon · sin(latMid) · sinc(dLat / 2), as `area`. + * Since the dLon of a closed ring sum to 0, sin(latMid) · sinc is replaced by its difference + * with sin(lat0), computed without cancellation, so that the error scales with the size of + * the ring, not with |sin(lat)|. + */ +function localArea(c: Float64Array, n: number, ctx: AreaContext) : number { + const { s0, c0, phi0 } = ctx; + let sum = 0; + let px = c[2 * n - 2], py = c[2 * n - 1]; + for (let i = 0; i < 2 * n; i += 2) { + const x = c[i], y = c[i + 1]; + const dl = x - px; + if (dl !== 0) { + // u = (latMid - lat0) / 2 in radians; sin(lat0 + 2u) - sin(lat0) = 2 cos(lat0 + u) sin(u) + const u = (y + py) * (HALF_DEGREES_TO_RADIANS / 2); + let su: number, cosLat: number; + if (u < 0.02 && u > -0.02) { + const u2 = u * u; + su = u * (1 - u2 / 6 * (1 - u2 / 20 * (1 - u2 / 42 * (1 - u2 / 72)))); + const cu = 1 - u2 / 2 * (1 - u2 / 12 * (1 - u2 / 30 * (1 - u2 / 56))); + cosLat = c0 * cu - s0 * su; + } else { + su = Math.sin(u); + cosLat = Math.cos(phi0 + u); + } + const dS = 2 * cosLat * su; + // sinc(h) - 1, h = dLat / 2 + const h = (y - py) * HALF_DEGREES_TO_RADIANS; + const h2 = h * h; + const sincm1 = h2 < 1e-4 + ? -h2 / 6 * (1 - h2 / 20 * (1 - h2 / 42 * (1 - h2 / 72))) + : Math.sin(h) / h - 1; + sum += dl * (dS * (1 + sincm1) + s0 * sincm1); + } + px = x; + py = y; + } + return -sum * DEGREES_TO_RADIANS; +} + +/** Unit-sphere area of a box in local coordinates. */ +function boxArea(x0: number, y0: number, x1: number, y1: number, ctx: AreaContext) : number { + const mid = ctx.phi0 + (y0 + y1) * HALF_DEGREES_TO_RADIANS; + return (x1 - x0) * DEGREES_TO_RADIANS * 2 * Math.cos(mid) * Math.sin((y1 - y0) * HALF_DEGREES_TO_RADIANS); +} + +function sumArea(rings: Ring[], ctx: AreaContext) : number { + let total = 0; + for (const r of rings) total += r.mu * localArea(r.c, r.n, ctx); + return total; +} + +// --- Axis-aligned Sutherland-Hodgman --- + +/** + * Keep the part of ring `src` (n positions) on one side of the line coordinate[axis] = v: + * `keepLow` keeps <= v, otherwise >= v. The crossing point is interpolated from the + * endpoint with the smaller coordinate, so that an edge shared by two rings gives the + * same point whatever its direction. Returns the number of positions written to `dst`. + */ +function axisPass(src: Float64Array, n: number, dst: Float64Array, axis: Axis, v: number, keepLow: boolean) : number { + const other = 1 - axis; + let m = 0; + let pa = src[2 * n - 2 + axis], pb = src[2 * n - 2 + other]; + for (let i = 0; i < 2 * n; i += 2) { + const qa = src[i + axis], qb = src[i + other]; + if ((pa < v && qa > v) || (pa > v && qa < v)) { + const xb = pa < qa ? pb + (v - pa) / (qa - pa) * (qb - pb) : qb + (v - qa) / (pa - qa) * (pb - qb); + if (axis === 0) { + dst[m++] = v; + dst[m++] = xb; + } else { + dst[m++] = xb; + dst[m++] = v; + } + } + if (keepLow ? qa <= v : qa >= v) { + dst[m++] = src[i]; + dst[m++] = src[i + 1]; + } + pa = qa; + pb = qb; + } + return m >> 1; +} + +/** Clip a ring to a box, or null when fewer than 3 positions survive. Only the sides crossing the ring's bbox are clipped. */ +function clipRingToBox(r: Ring, bx0: number, by0: number, bx1: number, by1: number) : Ring | null { + let src = r.c, n = r.n; + let toA = true; // output buffer: bufferA, then bufferB, alternately + const pass = (axis: Axis, v: number, keepLow: boolean) => { + const dst = toA ? growA(4 * n) : growB(4 * n); + toA = !toA; + n = axisPass(src, n, dst, axis, v, keepLow); + src = dst; + return n >= 3; + }; + if (r.x0 < bx0 && !pass(0, bx0, false)) return null; + if (r.x1 > bx1 && !pass(0, bx1, true)) return null; + if (r.y0 < by0 && !pass(1, by0, false)) return null; + if (r.y1 > by1 && !pass(1, by1, true)) return null; + return makeRing(src, 2 * n, r.mu, r.o); +} + +/** Split a ring by the line coordinate[axis] = v into its closed halves. */ +function splitRing(r: Ring, axis: Axis, v: number, low: Ring[], high: Ring[]) { + const c = r.c, n = r.n; + const lowBuffer = growA(4 * n), highBuffer = growB(4 * n); + const other = 1 - axis; + let nl = 0, nh = 0; + let pa = c[2 * n - 2 + axis], pb = c[2 * n - 2 + other]; + for (let i = 0; i < 2 * n; i += 2) { + const qa = c[i + axis], qb = c[i + other]; + if ((pa < v && qa > v) || (pa > v && qa < v)) { + const xb = pa < qa ? pb + (v - pa) / (qa - pa) * (qb - pb) : qb + (v - qa) / (pa - qa) * (pb - qb); + const x = axis === 0 ? v : xb, y = axis === 0 ? xb : v; + lowBuffer[nl++] = x; + lowBuffer[nl++] = y; + highBuffer[nh++] = x; + highBuffer[nh++] = y; + } + if (qa <= v) { + lowBuffer[nl++] = c[i]; + lowBuffer[nl++] = c[i + 1]; + } + if (qa >= v) { + highBuffer[nh++] = c[i]; + highBuffer[nh++] = c[i + 1]; + } + pa = qa; + pb = qb; + } + if (nl >= 6) low.push(makeRing(lowBuffer, nl, r.mu, r.o)); + if (nh >= 6) high.push(makeRing(highBuffer, nh, r.mu, r.o)); +} + +function splitRings(rings: Ring[], axis: Axis, v: number, low: Ring[], high: Ring[]) { + for (const r of rings) { + const lo = axis === 0 ? r.x0 : r.y0; + const hi = axis === 0 ? r.x1 : r.y1; + if (hi <= v) low.push(r); + else if (lo >= v) high.push(r); + else splitRing(r, axis, v, low, high); + } +} + +/** Rings clipped to a box; a ring inside the box is kept as is (rings are never mutated). */ +function clipRings(rings: Ring[], x0: number, y0: number, x1: number, y1: number) : Ring[] { + const out: Ring[] = []; + for (const r of rings) { + if (r.x1 < x0 || r.x0 > x1 || r.y1 < y0 || r.y0 > y1) continue; + if (r.x0 >= x0 && r.x1 <= x1 && r.y0 >= y0 && r.y1 <= y1) { + out.push(r); + continue; + } + const piece = clipRingToBox(r, x0, y0, x1, y1); + if (piece) out.push(piece); + } + return out; +} + +// --- Constant multiplicities --- + +/** + * Winding number inside the box of a ring lying on the box boundary (every edge along one + * side, compared exactly), or NaN when the ring does not lie on the boundary. + */ +function boundaryWinding(r: Ring, bx0: number, by0: number, bx1: number, by1: number) : number { + const c = r.c, n2 = 2 * r.n; + let px = c[n2 - 2], py = c[n2 - 1]; + for (let i = 0; i < n2; i += 2) { + const x = c[i], y = c[i + 1]; + if (!((x === bx0 && px === bx0) || (x === bx1 && px === bx1) || (y === by0 && py === by0) || (y === by1 && py === by1))) { + return NaN; + } + px = x; + py = y; + } + return Math.round(planarArea2(c, r.n) / (2 * (bx1 - bx0) * (by1 - by0))); +} + +/** Move the rings lying on the box boundary into a constant multiplicity, returned; the others go to `rest`. */ +function extractConstant(rings: Ring[], bx0: number, by0: number, bx1: number, by1: number, rest: Ring[]) : number { + let constant = 0; + for (const r of rings) { + const k = boundaryWinding(r, bx0, by0, bx1, by1); + if (Number.isNaN(k)) rest.push(r); + else constant += r.mu * k; + } + return constant; +} + +/** Multiplicity inside the box when every ring lies on its boundary, else NaN. */ +function constantMultiplicity(rings: Ring[], bx0: number, by0: number, bx1: number, by1: number) : number { + let constant = 0; + for (const r of rings) { + const k = boundaryWinding(r, bx0, by0, bx1, by1); + if (Number.isNaN(k)) return NaN; + constant += r.mu * k; + } + return constant; +} + +// --- Recursive split of the box --- + +/** Integral over the box of mF · mG, where mF (resp. mG) is sum(mu · winding) over the pieces of F (resp. G), all inside the box. */ +function boxIntegral(F: Ring[], G: Ring[], x0: number, y0: number, x1: number, y1: number, depth: number, ctx: AreaContext) : number { + if (F.length === 0 || G.length === 0) return 0; + const F2: Ring[] = [], G2: Ring[] = []; + const cF = extractConstant(F, x0, y0, x1, y1, F2); + const cG = extractConstant(G, x0, y0, x1, y1, G2); + let total = 0; + if (cF !== 0 && cG !== 0) total += cF * cG * boxArea(x0, y0, x1, y1, ctx); + if (cF !== 0 && G2.length) total += cF * sumArea(G2, ctx); + if (cG !== 0 && F2.length) total += cG * sumArea(F2, ctx); + if (F2.length && G2.length) total += solve(F2, G2, x0, y0, x1, y1, depth, ctx); + return total; +} + +/** Integral over the box of mF · mG, for non-empty sets of pieces not lying on the boundary. */ +function solve(F: Ring[], G: Ring[], x0: number, y0: number, x1: number, y1: number, depth: number, ctx: AreaContext) : number { + const nF = countPositions(F), nG = countPositions(G); + const small = Math.min(nF, nG), big = Math.max(nF, nG); + if (depth >= MAX_DEPTH || small * big <= Math.max(LEAF_PRODUCT, LEAF_FACTOR * (small + big))) { + return leaf(F, G, nF, nG, ctx); + } + const FL: Ring[] = [], FH: Ring[] = [], GL: Ring[] = [], GH: Ring[] = []; + if (x1 - x0 >= y1 - y0) { + const v = (x0 + x1) / 2; + if (!(v > x0 && v < x1)) return leaf(F, G, nF, nG, ctx); // box too small to split + splitRings(F, 0, v, FL, FH); + splitRings(G, 0, v, GL, GH); + return boxIntegral(FL, GL, x0, y0, v, y1, depth + 1, ctx) + boxIntegral(FH, GH, v, y0, x1, y1, depth + 1, ctx); + } + const v = (y0 + y1) / 2; + if (!(v > y0 && v < y1)) return leaf(F, G, nF, nG, ctx); + splitRings(F, 1, v, FL, FH); + splitRings(G, 1, v, GL, GH); + return boxIntegral(FL, GL, x0, y0, x1, v, depth + 1, ctx) + boxIntegral(FH, GH, x0, v, x1, y1, depth + 1, ctx); +} + +// --- Leaves: triangulate the smaller side, clip the other side's rings by each triangle --- + +function leaf(F: Ring[], G: Ring[], nF: number, nG: number, ctx: AreaContext) : number { + const [A, B, nA, nB]: [Ring[], Ring[], number, number] = nF <= nG ? [F, G, nF, nG] : [G, F, nG, nF]; + // About one triangle per position of A, each clipping every position of B. + if (nA * nB > MAX_LEAF_OPERATIONS) throw FALLBACK; + let total = 0; + for (const a of A) total += a.mu * ringIntegral(a, B, ctx); + return total; +} + +/** Integral of winding(a) · mB, where mB = sum(mu · winding) over the rings of B. */ +function ringIntegral(a: Ring, B: Ring[], ctx: AreaContext) : number { + const c = a.c, n = a.n; + const area2 = planarArea2(c, n); + if (area2 === 0) return 0; + if (n === 3) return triangleIntegral(B, c[0], c[1], c[2], c[3], c[4], c[5], ctx) * (area2 > 0 ? 1 : -1); + const indices = earcut(c, null, 2); + if (deviation(c, null, 2, indices) <= MAX_DEVIATION) { + // earcut triangles cover the inside of `a` once, where winding(a) = sign(area). + let total = 0; + for (let i = 0; i < indices.length; i += 3) { + const i0 = 2 * indices[i], i1 = 2 * indices[i + 1], i2 = 2 * indices[i + 2]; + total += triangleIntegral(B, c[i0], c[i0 + 1], c[i1], c[i1 + 1], c[i2], c[i2 + 1], ctx); + } + return area2 > 0 ? total : -total; + } + // earcut fails on the zero-width bridges left by Sutherland-Hodgman, which a signed fan + // handles exactly, but also on self-intersecting rings, left to `@turf/intersect`. + if (hasProperCrossing(c, n)) throw FALLBACK; + // Signed fan from the first position: winding(a) = sum of the signed windings of its triangles. + let total = 0; + const ax = c[0], ay = c[1]; + for (let i = 2; i + 3 < 2 * n; i += 2) { + const bx = c[i], by = c[i + 1], cx = c[i + 2], cy = c[i + 3]; + const cr = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax); + if (cr > 0) total += triangleIntegral(B, ax, ay, bx, by, cx, cy, ctx); + else if (cr < 0) total -= triangleIntegral(B, ax, ay, bx, by, cx, cy, ctx); + } + return total; +} + +/** + * True when two edges of the ring cross at a point interior to both. Zero-width bridges + * (collinear overlaps along a clipping line) and touching positions are not crossings, so + * the pieces of a simple ring have none. + */ +function hasProperCrossing(c: Float64Array, n: number) : boolean { + for (let i = 0; i < n; i++) { + const i2 = 2 * i, j2 = i + 1 < n ? i2 + 2 : 0; + const ax = c[i2], ay = c[i2 + 1], bx = c[j2], by = c[j2 + 1]; + const minX = Math.min(ax, bx), maxX = Math.max(ax, bx), minY = Math.min(ay, by), maxY = Math.max(ay, by); + for (let k = i + 2; k < n; k++) { + if (i === 0 && k === n - 1) continue; // adjacent through the closing edge + const k2 = 2 * k, l2 = k + 1 < n ? k2 + 2 : 0; + const cx = c[k2], cy = c[k2 + 1], dx = c[l2], dy = c[l2 + 1]; + if (Math.max(cx, dx) < minX || Math.min(cx, dx) > maxX || Math.max(cy, dy) < minY || Math.min(cy, dy) > maxY) continue; + const o1 = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax); + const o2 = (bx - ax) * (dy - ay) - (by - ay) * (dx - ax); + if (!((o1 > 0 && o2 < 0) || (o1 < 0 && o2 > 0))) continue; + const o3 = (dx - cx) * (ay - cy) - (dy - cy) * (ax - cx); + const o4 = (dx - cx) * (by - cy) - (dy - cy) * (bx - cx); + if ((o3 > 0 && o4 < 0) || (o3 < 0 && o4 > 0)) return true; + } + } + return false; +} + +/** Integral of mB over a triangle (any orientation), i.e. sum(mu · signed area of ring ∩ triangle). */ +function triangleIntegral(B: Ring[], ax: number, ay: number, bx: number, by: number, cx: number, cy: number, ctx: AreaContext) : number { + const cr = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax); + if (cr === 0) return 0; + if (cr < 0) { + // Counterclockwise, so that the interior lies on the left of each edge. + const tx = bx, ty = by; + bx = cx; + by = cy; + cx = tx; + cy = ty; + } + const tx0 = Math.min(ax, bx, cx), tx1 = Math.max(ax, bx, cx); + const ty0 = Math.min(ay, by, cy), ty1 = Math.max(ay, by, cy); + // A piece of a simple ring has the orientation of the ring: an opposite one, well beyond + // rounding, reveals a self-intersecting ring. The rounding of `localArea` is a few ulps of + // |dLon| · |local latitude| per edge, hence the second term. + const extent = Math.max(tx1 - tx0, ty1 - ty0); + const magnitude = Math.max(extent, Math.abs(tx0), Math.abs(tx1), Math.abs(ty0), Math.abs(ty1)); + const relative = 1e-6 * (tx1 - tx0) * (ty1 - ty0) * DEGREES_TO_RADIANS * DEGREES_TO_RADIANS; + const rounding = 1e-12 * extent * magnitude * DEGREES_TO_RADIANS * DEGREES_TO_RADIANS; + let total = 0; + for (const r of B) { + if (r.x1 < tx0 || r.x0 > tx1 || r.y1 < ty0 || r.y0 > ty1) continue; + const clipped = clippedArea(r, ax, ay, bx, by, cx, cy, ctx); + if (clipped * r.o < -(relative + rounding * r.n)) throw FALLBACK; + total += r.mu * clipped; + } + return total; +} + +/** + * Sutherland-Hodgman pass keeping the left of the line s -> e. The crossing point is + * interpolated from the endpoint lying strictly inside, whatever the edge direction. + */ +function trianglePass(src: Float64Array, n: number, dst: Float64Array, sx: number, sy: number, ex: number, ey: number) : number { + const dx = ex - sx, dy = ey - sy; + let m = 0; + let px = src[2 * n - 2], py = src[2 * n - 1]; + let ps = dx * (py - sy) - dy * (px - sx); + for (let i = 0; i < 2 * n; i += 2) { + const qx = src[i], qy = src[i + 1]; + const qs = dx * (qy - sy) - dy * (qx - sx); + if (ps > 0 ? qs < 0 : ps < 0 && qs > 0) { + if (ps > 0) { + const t = ps / (ps - qs); + dst[m++] = px + t * (qx - px); + dst[m++] = py + t * (qy - py); + } else { + const t = qs / (qs - ps); + dst[m++] = qx + t * (px - qx); + dst[m++] = qy + t * (py - qy); + } + } + if (qs >= 0) { + dst[m++] = qx; + dst[m++] = qy; + } + px = qx; + py = qy; + ps = qs; + } + return m >> 1; +} + +/** Unit-sphere signed area of ring ∩ triangle (counterclockwise triangle). */ +function clippedArea(r: Ring, ax: number, ay: number, bx: number, by: number, cx: number, cy: number, ctx: AreaContext) : number { + let n = r.n; + let dst = growA(4 * n); + n = trianglePass(r.c, n, dst, ax, ay, bx, by); + if (n < 3) return 0; + let src = dst; + dst = growB(4 * n); + n = trianglePass(src, n, dst, bx, by, cx, cy); + if (n < 3) return 0; + src = dst; + dst = growA(4 * n); + n = trianglePass(src, n, dst, cx, cy, ax, ay); + if (n < 3) return 0; + return localArea(dst, n, ctx); +} + +// --- Prepared filter --- + +type PreparedFilter = { + /** Local origin: the centre of the filter bbox. */ + ox: number; + oy: number; + /** True when a filter ring has no area but some extent (self-intersecting, like a bow tie). */ + invalid: boolean; + /** Rings of the tile containing an absolute bbox, in local coordinates. */ + near: (minX: number, minY: number, maxX: number, maxY: number) => Ring[]; +}; + +/** + * Flat copy of the filter rings, clipped lazily to nested tiles, each from its parent. Tiles + * are 1.5 times as wide as their cell, so that a bbox at most half a cell wide lies in the + * tile of the cell containing its south-west corner; a ring lying inside a tile is shared, + * not copied. + */ +function prepareFilter(filter: Polygon | MultiPolygon) : PreparedFilter { + const [west, south, east, north] = bbox(filter); + const ox = (west + east) / 2, oy = (south + north) / 2; + const size = Math.max(east - west, north - south); + const root: Ring[] = []; + let invalid = false; + for (const polygon of polygonsOf(filter)) { + for (let index = 0; index < polygon.length; index++) { + const ring = localRing(polygon[index], ox, oy); + if (!ring) continue; + const area2 = planarArea2(ring.c, ring.n); + if (Math.abs(area2) <= 1e-12 * 2 * (ring.x1 - ring.x0) * (ring.y1 - ring.y0)) { + if (ring.x1 > ring.x0 && ring.y1 > ring.y0) invalid = true; + continue; + } + ring.o = area2 > 0 ? 1 : -1; + ring.mu = (index == 0 ? 1 : -1) * ring.o; + root.push(ring); + } + } + const localWest = west - ox, localSouth = south - oy; + const tiles = new Map(); + function tile(zoom: number, x: number, y: number) : Ring[] { + if (zoom == 0) return root; + const key = `${zoom}/${x}/${y}`; + let rings = tiles.get(key); + if (rings === undefined) { + const parent = tile(zoom - 1, Math.floor(x / 2), Math.floor(y / 2)); + const cell = size / 2 ** zoom; + rings = parent.length == 0 ? parent : clipRings(parent, localWest + x * cell, localSouth + y * cell, localWest + (x + 1.5) * cell, localSouth + (y + 1.5) * cell); + tiles.set(key, rings); + } + return rings; + } + return { + ox, oy, invalid, + near(minX, minY, maxX, maxY) { + // Deepest tile whose cell is at least twice as wide as the bbox. + const width = Math.max(maxX - minX, maxY - minY); + let zoom = 0; + while (zoom < 20 && 2 * width <= size / 2 ** (zoom + 1)) zoom++; + const cell = size / 2 ** zoom; + return tile(zoom, Math.floor((minX - west) / cell), Math.floor((minY - south) / cell)); + }, + }; +} + +// --- Entry point --- + +/** + * Prepare the filter once (per request), and return the function computing the area (m²) + * of the part of each feature lying inside it. + */ +export function makeIntersectionArea(filter: Polygon | MultiPolygon) : IntersectionArea { + const prepared = prepareFilter(filter); + return (geo) => { + const featureBbox = bbox(geo); + if (prepared.invalid) return filterFallback(filter, geo, featureBbox); + // The filter near the feature: its tile, then exactly the feature bbox. + const [f0, f1, f2, f3] = featureBbox; + const { ox, oy } = prepared; + const fx0 = f0 - ox, fy0 = f1 - oy, fx1 = f2 - ox, fy1 = f3 - oy; + let near: Ring[]; + try { + near = clipRings(prepared.near(f0, f1, f2, f3), fx0, fy0, fx1, fy1); + } catch (error) { + // A tile piece with the wrong orientation: self-intersecting filter ring. + if (error !== FALLBACK) throw error; + return filterFallback(filter, geo, featureBbox); + } + if (near.length == 0) return 0; + let W0 = Infinity, W1 = Infinity, W2 = -Infinity, W3 = -Infinity; + for (const r of near) { + if (r.x0 < W0) W0 = r.x0; + if (r.y0 < W1) W1 = r.y0; + if (r.x1 > W2) W2 = r.x1; + if (r.y1 > W3) W3 = r.y1; + } + if (!(W2 > W0 && W3 > W1)) return 0; // no areal overlap + if (W0 === fx0 && W1 === fy0 && W2 === fx1 && W3 === fy1) { + // The filter covers the whole feature bbox, or none of it. + const constant = constantMultiplicity(near, W0, W1, W2, W3); + const result = Number.isNaN(constant) ? null : insideShortcut(constant, geo, featureBbox); + if (result !== null) return result; + } + const fallbackInputs = () : [Polygon | MultiPolygon, BBox] | null => { + const clippedFilter = dropEmptyRings(bboxClip(filter, featureBbox).geometry); + return clippedFilter ? [clippedFilter, bbox(clippedFilter)] : null; + }; + return featureIntegral(geo, near, W0, W1, W2, W3, ox, oy, fallbackInputs); + }; +} + +/** `@turf/intersect` on the feature and the whole filter clipped to its bbox. */ +function filterFallback(filter: Polygon | MultiPolygon, geo: Polygon | MultiPolygon, featureBbox: BBox) : number { + const clippedFilter = dropEmptyRings(bboxClip(filter, featureBbox).geometry); + if (!clippedFilter) return 0; + return snap(polyclipFallback(geo, clippedFilter, bbox(clippedFilter)), geo, area(geo), false); +} + +/** `@turf/intersect` on the pre-clipped pair, one feature polygon at a time, so that polygons are summed. */ +function polyclipFallback(geo: Polygon | MultiPolygon, clippedFilter: Polygon | MultiPolygon, workingBox: BBox) : number { + let total = 0; + for (const polygon of polygonsOf(geo)) { + const clipped = dropEmptyRings(bboxClip({ type: "Polygon", coordinates: polygon }, workingBox).geometry); + if (!clipped) continue; + const inter = intersect(featureCollection([feature(clipped), feature(clippedFilter)])); + total += inter == null ? 0 : area(inter.geometry); + } + return total; +} + +/** + * Result when the filter has a constant multiplicity over the whole feature bbox, or null to + * take the general path: a feature without area but with some extent may be a bow tie, + * whose `area` is not the area it covers. + */ +function insideShortcut(constant: number, geo: Polygon | MultiPolygon, [f0, f1, f2, f3]: BBox) : number | null { + if (constant <= 0) return 0; + const featureArea = area(geo); + const bboxArea = (f2 - f0) * (f3 - f1) * DEGREES_TO_RADIANS * DEGREES_TO_RADIANS * Math.cos((f1 + f3) / 2 * DEGREES_TO_RADIANS) * R2; + if (!(featureArea > 1e-12 * bboxArea)) return null; + return featureArea; +} + +/** + * Intersection area of the feature with the filter rings G, all in local coordinates + * (relative to lon0, lat0), G lying inside the working box. + */ +function featureIntegral( + geo: Polygon | MultiPolygon, + G: Ring[], + bx0: number, by0: number, bx1: number, by1: number, + lon0: number, lat0: number, + fallbackInputs: () => [Polygon | MultiPolygon, BBox] | null, +) : number { + const phi0 = lat0 * DEGREES_TO_RADIANS; + const ctx: AreaContext = { phi0, s0: Math.sin(phi0), c0: Math.cos(phi0) }; + + // The feature clipped to the working box. + const F: Ring[] = []; + const whole: Ring[] = []; // every ring of the feature, unclipped, for the snapping decision + let clipped = false; + let invalid = false; + let planarFeature2 = 0; // twice the planar area of the feature, for the snapping estimate + for (const polygon of polygonsOf(geo)) { + for (let index = 0; index < polygon.length; index++) { + const ring = localRing(polygon[index], lon0, lat0); + if (!ring) continue; + const area2 = planarArea2(ring.c, ring.n); + const bboxArea2 = 2 * (ring.x1 - ring.x0) * (ring.y1 - ring.y0); + if (Math.abs(area2) <= 1e-12 * bboxArea2) { + // No area but some extent: a self-intersecting ring whose lobes cancel (bow tie) + // is left to `@turf/intersect`; a zero-width ring is only slower there. + if (bboxArea2 > 0) invalid = true; + continue; + } + const role = index == 0 ? 1 : -1; + ring.o = area2 > 0 ? 1 : -1; + ring.mu = role * ring.o; + planarFeature2 += role * Math.abs(area2); + whole.push(ring); + if (ring.x0 >= bx0 && ring.x1 <= bx1 && ring.y0 >= by0 && ring.y1 <= by1) { + F.push(ring); + continue; + } + clipped = true; + if (ring.x1 < bx0 || ring.x0 > bx1 || ring.y1 < by0 || ring.y0 > by1 || invalid) continue; + try { + const piece = clipRingToBox(ring, bx0, by0, bx1, by1); + if (piece) F.push(piece); + } catch (error) { + if (error !== FALLBACK) throw error; + invalid = true; // a piece with the wrong orientation: self-intersecting ring + } + } + } + const estimate = planarFeature2 / 2 * DEGREES_TO_RADIANS * DEGREES_TO_RADIANS * ctx.c0 * R2; + + if (!invalid) { + try { + const F2: Ring[] = [], G2: Ring[] = []; + const cF = extractConstant(F, bx0, by0, bx1, by1, F2); + const cG = extractConstant(G, bx0, by0, bx1, by1, G2); + if (G2.length == 0 && !clipped) { + // The filter covers the whole working box, which contains the feature, or none of it. + return cG > 0 ? area(geo) : 0; + } + let total = 0; + if (cF !== 0 && cG !== 0) total += cF * cG * boxArea(bx0, by0, bx1, by1, ctx); + if (cF !== 0 && G2.length) total += cF * sumArea(G2, ctx); + if (cG !== 0 && F2.length) total += cG * sumArea(F2, ctx); + if (F2.length && G2.length) total += solve(F2, G2, bx0, by0, bx1, by1, 0, ctx); + // Capped at the feature area. + return snap(total * R2, geo, estimate, true, () => sumArea(whole, ctx) * R2); + } catch (error) { + if (error !== FALLBACK) throw error; + } + } + // Not capped: for a self-intersecting feature, `area` is not the area it covers. + const inputs = fallbackInputs(); + if (!inputs) return 0; + return snap(polyclipFallback(geo, inputs[0], inputs[1]), geo, estimate, false); +} + +/** + * Snap to the feature area (within SNAP_FULL, or above it when `cap`) or to 0 (within + * SNAP_ZERO). The exact feature area is only computed when the planar estimate says that + * the result may be close to either. + */ +function snap(total: number, geo: Polygon | MultiPolygon, estimate: number, cap: boolean, localFeatureArea?: () => number) : number { + if (Number.isNaN(total)) throw new Error("intersection area is NaN"); + if (total <= 0) return 0; + if (total <= 100 * SNAP_ZERO * estimate || total >= 0.5 * estimate) { + // Decided on the same (local) formula as `total` when available; the snapped value is + // `area`, so that a feature inside the filter gets exactly its `area`. + const featureArea = localFeatureArea ? localFeatureArea() : area(geo); + if (total >= featureArea * (1 - SNAP_FULL) && (cap || total <= featureArea * (1 + SNAP_FULL))) return area(geo); + if (total <= featureArea * SNAP_ZERO) return 0; + } + return total; +} diff --git a/src/wfs/spatialExtras.ts b/src/wfs/spatialExtras.ts index 0fec7650..0bc8ebc7 100644 --- a/src/wfs/spatialExtras.ts +++ b/src/wfs/spatialExtras.ts @@ -1,12 +1,11 @@ import { centroid } from "@turf/centroid"; import { bbox } from "@turf/bbox"; import turfLength from "@turf/length"; -import { intersect } from "@turf/intersect"; import { circle } from "@turf/circle"; import { bboxPolygon } from "@turf/bbox-polygon"; -import type { BBox, Geometry, LineString, MultiLineString, MultiPolygon, Point, Polygon, Position } from "geojson"; +import type { Geometry, LineString, MultiLineString, MultiPolygon, Point, Polygon, Position } from "geojson"; import distance from "../helpers/distance.js"; -import { feature, featureCollection } from "@turf/helpers"; +import { feature } from "@turf/helpers"; import { getSpatialFilter } from "./spatialFilter.js"; import type { SpatialFilterInput, @@ -14,7 +13,9 @@ import type { SpatialFilter, } from "./schema.js"; import { bboxClip } from "@turf/bbox-clip"; -import area, { clipRingToTriangle, triangulate } from "../helpers/area.js"; +import area from "../helpers/area.js"; +import { dropEmptyRings } from "../helpers/geojson.js"; +import { makeIntersectionArea, type IntersectionArea } from "../helpers/intersectionArea.js"; export type FeatureCollectionPostProcessInput = { typename: string, @@ -134,7 +135,7 @@ function isComputableGeometry(geometry: unknown) : geometry is Geometry { type SpatialContext = { filterCentroid: Point | null, - clipFilter: FilterClipper | null + intersectionArea: IntersectionArea | null } export function prepareSpatialContext(input: FeatureCollectionPostProcessInput, resolvedGeometryRef?: Geometry) : SpatialContext { @@ -143,7 +144,7 @@ export function prepareSpatialContext(input: FeatureCollectionPostProcessInput, const context : SpatialContext = { filterCentroid: null, - clipFilter: null, + intersectionArea: null, }; if (!requires_distance_to_filter_center && !requires_intersection_area) { @@ -165,161 +166,20 @@ export function prepareSpatialContext(input: FeatureCollectionPostProcessInput, if (requires_intersection_area) { try { const filterPolygons = geometryToPolygons(spatialFilterToGeometry(spatialFilter, resolvedGeometryRef)); - context.clipFilter = filterPolygons && makeFilterClipper(filterPolygons); + context.intersectionArea = filterPolygons && makeIntersectionArea(filterPolygons); } catch {} } return context; } -/** 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; -} - -/** Return the filter clipped to a tile containing the bbox, or null if nothing areal remains. */ -type FilterClipper = (bbox: BBox) => Polygon | MultiPolygon | null; - -/** Build a `FilterClipper` which only processes the filter vertices near the bbox. - * - * The filter is clipped to nested tiles, each from the cached clip of its parent. - * Tiles are 1.5 times as wide as their cell, so that a bbox at most half a cell - * wide lies in the tile of the cell containing its south-west corner. - */ -function makeFilterClipper(filter: Polygon | MultiPolygon) : FilterClipper { - const [west, south, east, north] = bbox(filter); - const size = Math.max(east - west, north - south); - const tiles = new Map(); - function tile(zoom: number, x: number, y: number) : Polygon | MultiPolygon | null { - if (zoom == 0) return filter; - const key = `${zoom}/${x}/${y}`; - if (!tiles.has(key)) { - const parent = tile(zoom - 1, Math.floor(x / 2), Math.floor(y / 2)); - const cell = size / 2 ** zoom; - tiles.set(key, parent && dropEmptyRings(bboxClip(parent, [west + x * cell, south + y * cell, west + (x + 1.5) * cell, south + (y + 1.5) * cell]).geometry)); - } - return tiles.get(key) ?? null; - } - return ([minX, minY, maxX, maxY]) => { - // Deepest tile whose cell is at least twice as wide as the bbox. - const width = Math.max(maxX - minX, maxY - minY); - let zoom = 0; - while (zoom < 20 && 2 * width <= size / 2 ** (zoom + 1)) zoom++; - const cell = size / 2 ** zoom; - return tile(zoom, Math.floor((minX - west) / cell), Math.floor((minY - south) / cell)); - }; -} - -/** Polygons of a polygonal geometry, as lists of rings. */ -function polygonsOf(geo: Polygon | MultiPolygon) : Position[][][] { - return geo.type == "Polygon" ? [geo.coordinates] : geo.coordinates; -} - -function positionCount(geo: Polygon | MultiPolygon) : number { - return polygonsOf(geo).flat().reduce((count, ring) => count + ring.length, 0); -} - -/** True when `geo` is the whole `box`, as clipping returns a geometry covering it. */ -function isBox(geo: Polygon | MultiPolygon, [west, south, east, north]: BBox) : boolean { - const rings = polygonsOf(geo).flat(); - if (rings.length != 1 || rings[0].length != 5) return false; - const corners = rings[0].slice(0, -1); - return corners.every(([x, y]) => (x == west || x == east) && (y == south || y == north)) - && new Set(corners.map(String)).size == 4; -} - -/** Beyond this many positions, geometries are split before being triangulated. */ -const MAX_TRIANGULATED_POSITIONS = 32; -/** Bound on the number of halvings, in case positions pile up at the same place. */ -const MAX_HALVINGS = 24; - -/** Return the area of the intersection between two areal geometries lying in `box`. - * - * Clipping by a triangle is simple, fast and robust, unlike general polygon - * clipping: the geometry with fewer positions is split into triangles, by which - * the other one is clipped. As this costs the product of their numbers of - * positions, the box is first halved until one of them is small. - */ -function boxIntersectionArea(a: Polygon | MultiPolygon, b: Polygon | MultiPolygon, box: BBox, halvings = 0) : number { - const [small, large] = positionCount(a) <= positionCount(b) ? [a, b] : [b, a]; - if (positionCount(small) > MAX_TRIANGULATED_POSITIONS && halvings < MAX_HALVINGS) { - const [west, south, east, north] = box; - const halves : BBox[] = east - west > north - south - ? [[west, south, (west + east) / 2, north], [(west + east) / 2, south, east, north]] - : [[west, south, east, (south + north) / 2], [west, (south + north) / 2, east, north]]; - return halves.reduce((total, half) => { - const smallHalf = dropEmptyRings(bboxClip(small, half).geometry); - const largeHalf = smallHalf && dropEmptyRings(bboxClip(large, half).geometry); - return largeHalf ? total + boxIntersectionArea(smallHalf, largeHalf, half, halvings + 1) : total; - }, 0); - } - - if (isBox(small, box)) return area(large); - const triangulations = polygonsOf(small).map(triangulate); - if (triangulations.includes(null)) { - // Fall back on general polygon clipping, as for self-intersecting rings. - const inter = intersect(featureCollection([feature(small), feature(large)])); - return inter == null ? 0 : area(inter.geometry); - } - let total = 0; - for (const triangle of triangulations.flatMap((triangles) => triangles!)) { - total += area({ type: "MultiPolygon", coordinates: polygonsOf(large).map((rings) => rings.map((ring) => clipRingToTriangle(ring, triangle))) }); - } - return total; -} - -/** Return the area of the intersection between two areal (2D) geometries. */ -function polygonsIntersectionArea(geo: Polygon | MultiPolygon, clipFilter: FilterClipper) : number { - let total = 0; - for (const polygon of polygonsOf(geo)) { - const polygonGeometry : Polygon = { type: "Polygon", coordinates: polygon }; - // Whatever lies outside the bbox of one geometry cannot intersect it, so - // clipping the other first leaves the result unchanged while only the - // neighbouring positions are processed afterwards. - const polygonBbox = bbox(polygonGeometry); - const nearFilter = clipFilter(polygonBbox); - const clippedFilter = nearFilter && dropEmptyRings(bboxClip(nearFilter, polygonBbox).geometry); - if (!clippedFilter) continue; // no areal overlap - const polygonArea = area(polygonGeometry); - if (isBox(clippedFilter, polygonBbox)) { // the filter contains the polygon - total += polygonArea; - continue; - } - const filterBbox = bbox(clippedFilter); - // Clipping the polygon to its own bbox would only copy it. - const clippedPolygon = filterBbox.every((value, index) => value == polygonBbox[index]) - ? polygonGeometry - : dropEmptyRings(bboxClip(polygonGeometry, filterBbox).geometry); - if (!clippedPolygon) continue; - - const covered = boxIntersectionArea(clippedPolygon, clippedFilter, filterBbox); - // Snap to 0 or to the whole polygon despite rounding errors. - total += covered > polygonArea * (1 - 1e-9) ? polygonArea : covered > polygonArea * 1e-9 ? covered : 0; - } - return total; -} - /** Return the area (m²) of the part of a geometry lying inside a spatial filter. * * null when it cannot be computed: the geometry or the filter has no areal part, * or the filter geometry could not be prepared. 0 only when both are areal and * do not overlap. */ -function intersectionAreaWithSpatialFilter(geom: Geometry, spatialFilter: SpatialFilter, clipFilter: FilterClipper | null) : number | null { +function intersectionAreaWithSpatialFilter(geom: Geometry, spatialFilter: SpatialFilter, intersectionArea: IntersectionArea | null) : number | null { const geo = geometryToPolygons(geom); if (!geo) { return null; // non-areal feature @@ -332,8 +192,8 @@ function intersectionAreaWithSpatialFilter(geom: Geometry, spatialFilter: Spatia case "dwithin_point": case "intersects_feature": case "travel_time": { - if (!clipFilter) return null; // non-areal filter, or filter preparation failed - return polygonsIntersectionArea(geo, clipFilter); + if (!intersectionArea) return null; // non-areal filter, or filter preparation failed + return intersectionArea(geo); } case "bbox": { const clipped = dropEmptyRings(bboxClip(geo, [spatialFilter.west, spatialFilter.south, spatialFilter.east, spatialFilter.north]).geometry); @@ -436,7 +296,7 @@ export function deriveFromGeometry(geometry: unknown, input: FeatureCollectionPo if (requires_intersection_area) { try { - ret.intersection_area = intersectionAreaWithSpatialFilter(geo, spatialFilter, context.clipFilter); + ret.intersection_area = intersectionAreaWithSpatialFilter(geo, spatialFilter, context.intersectionArea); } catch { ret.intersection_area = null; } diff --git a/test/helpers/area.test.ts b/test/helpers/area.test.ts index 55a02c4d..4c44af6b 100644 --- a/test/helpers/area.test.ts +++ b/test/helpers/area.test.ts @@ -1,17 +1,22 @@ import { describe, expect, it } from "vitest"; -import { clipRingToTriangle, sphericalRingArea } from "../../src/helpers/area.js"; +import { sphericalRingArea } from "../../src/helpers/area.js"; describe("helpers/area", () => { - describe("clipRingToTriangle", () => { - it("should clip by a triangle whatever its winding", () => { - // The triangle covers the north-eastern quarter of the square. - const square = [[0, 0], [2, 0], [2, 2], [0, 2], [0, 0]]; - const quarter = sphericalRingArea([[1, 1], [2, 1], [2, 2], [1, 2], [1, 1]]); - const [a, b, c] = [[1, 1], [3, 1], [1, 3]]; + describe("sphericalRingArea", () => { + it("should be exact for a lon/lat box", () => { + // dLon · (sin(north) - sin(south)) on the unit sphere + const box = [[2, 48], [3, 48], [3, 49], [2, 49], [2, 48]]; + const expected = (Math.PI / 180) * (Math.sin(49 * Math.PI / 180) - Math.sin(48 * Math.PI / 180)); - expect(sphericalRingArea(clipRingToTriangle(square, [a, b, c]))).toBeCloseTo(quarter, 12); - expect(sphericalRingArea(clipRingToTriangle(square, [a, c, b]))).toBeCloseTo(quarter, 12); + expect(sphericalRingArea(box)).toBeCloseTo(expected, 15); + }); + + it("should not change when an edge is split, as clipping does", () => { + const triangle = [[2, 48], [3, 48.2], [2.4, 49], [2, 48]]; + const split = [[2, 48], [2.5, 48.1], [3, 48.2], [2.4, 49], [2, 48]]; + + expect(sphericalRingArea(split)).toBeCloseTo(sphericalRingArea(triangle), 15); }); }); }); diff --git a/test/helpers/intersectionArea.test.ts b/test/helpers/intersectionArea.test.ts new file mode 100644 index 00000000..5579a202 --- /dev/null +++ b/test/helpers/intersectionArea.test.ts @@ -0,0 +1,69 @@ +import { describe, expect, it } from "vitest"; +import { intersect } from "@turf/intersect"; +import { feature, featureCollection } from "@turf/helpers"; +import type { MultiPolygon, Polygon, Position } from "geojson"; + +import area from "../../src/helpers/area.js"; +import { makeIntersectionArea } from "../../src/helpers/intersectionArea.js"; + +function rectangle(west: number, south: number, east: number, north: number) : Position[] { + return [[west, south], [east, south], [east, north], [west, north], [west, south]]; +} + +function wavyRing(centerLon: number, centerLat: number, radius: number, vertices: number, waves: number) : Position[] { + const positions = Array.from({ length: vertices }, (_, index) => { + const angle = (Math.PI * 2 * index) / vertices; + const wavyRadius = radius * (1 + 0.05 * Math.sin(waves * angle)); + return [centerLon + wavyRadius * Math.cos(angle), centerLat + wavyRadius * Math.sin(angle)]; + }); + return [...positions, positions[0]]; +} + +function polyclipArea(a: Polygon | MultiPolygon, b: Polygon | MultiPolygon) : number { + const inter = intersect(featureCollection([feature(a), feature(b)])); + return inter == null ? 0 : area(inter.geometry); +} + +describe("helpers/intersectionArea", () => { + const filter: Polygon = { type: "Polygon", coordinates: [rectangle(2, 48, 3, 49)] }; + + it("should sum the overlapping polygons of a MultiPolygon, as area does", () => { + const overlapping: MultiPolygon = { + type: "MultiPolygon", + coordinates: [[rectangle(2.2, 48.2, 2.6, 48.6)], [rectangle(2.4, 48.4, 2.8, 48.8)]], + }; + const intersectionArea = makeIntersectionArea(filter); + + expect(intersectionArea(overlapping)).toEqual(area(overlapping)); + + // Shifted across the eastern edge: each polygon is clipped on its own. + const crossing: MultiPolygon = { + type: "MultiPolygon", + coordinates: [[rectangle(2.8, 48.2, 3.2, 48.6)], [rectangle(2.9, 48.4, 3.3, 48.8)]], + }; + const expected = area({ type: "MultiPolygon", coordinates: [[rectangle(2.8, 48.2, 3, 48.6)], [rectangle(2.9, 48.4, 3, 48.8)]] }); + expect(intersectionArea(crossing)).toBeCloseTo(expected, 3); + }); + + it("should match general polygon clipping on detailed geometries crossing each other", () => { + const reference: Polygon = { type: "Polygon", coordinates: [wavyRing(2.35, 48.85, 0.3, 5000, 150)] }; + const intersectionArea = makeIntersectionArea(reference); + + for (const [lon, lat] of [[2.65, 48.85], [2.35, 48.55], [2.1, 49.05]]) { + const geo: Polygon = { type: "Polygon", coordinates: [wavyRing(lon, lat, 0.12, 3000, 90)] }; + const expected = polyclipArea(geo, reference); + + expect(Math.abs(intersectionArea(geo) - expected) / area(geo)).toBeLessThan(1e-9); + } + }); + + it("should fall back on general polygon clipping for a self-intersecting filter", () => { + const bowTie: Polygon = { + type: "Polygon", + coordinates: [[[2, 48], [3, 49], [3, 48], [2, 49], [2, 48]]], + }; + const geo: Polygon = { type: "Polygon", coordinates: [rectangle(2.1, 48.1, 2.9, 48.4)] }; + + expect(makeIntersectionArea(bowTie)(geo)).toBeCloseTo(polyclipArea(geo, bowTie), 0); + }); +}); diff --git a/test/wfs/geometry.test.ts b/test/wfs/geometry.test.ts index 3df380e9..b1fa34d2 100644 --- a/test/wfs/geometry.test.ts +++ b/test/wfs/geometry.test.ts @@ -2,7 +2,7 @@ import { describe, expect, it } from "vitest"; import { geometryToEwkt } from "../../src/wfs/geometry"; import type { Geometry } from "geojson"; import { isGeometryLike } from "../../src/helpers/geojson"; -import { dropEmptyRings } from "../../src/wfs/spatialExtras"; +import { dropEmptyRings } from "../../src/helpers/geojson"; describe("geometryToEwkt", () => { // --- Point and MultiPoint (already partially covered via queryPreparation tests) --- diff --git a/test/wfs/response.test.ts b/test/wfs/response.test.ts index cb2537c9..113fcd1b 100644 --- a/test/wfs/response.test.ts +++ b/test/wfs/response.test.ts @@ -639,7 +639,8 @@ describe("wfs_engine/response", () => { const inside = rectangle(2.4, 48.4, 2.6, 48.6); const feature = deriveExtras({ type: "MultiPolygon", coordinates: [[rectangle(3, 48.4, 3.2, 48.6)], [inside]] }); - expect(feature.intersection_area).toEqual(area({ type: "Polygon", coordinates: [inside] })); + // Not snapped to the area of the inside part: the MultiPolygon is integrated as a whole. + expect(feature.intersection_area).toBeCloseTo(area({ type: "Polygon", coordinates: [inside] }), 3); }); /** Area of the intersection computed by general polygon clipping, as a reference. */ @@ -669,6 +670,25 @@ describe("wfs_engine/response", () => { expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(withHole, referenceSquare), 0); }); + it("should fall back on general polygon clipping for a self-intersecting geometry", () => { + // Bow tie, whose lobes cancel in its signed area. + const bowTie: Polygon = { + type: "Polygon", + coordinates: [[[2.8, 48.4], [3.2, 48.6], [3.2, 48.4], [2.8, 48.6], [2.8, 48.4]]], + }; + const feature = deriveExtras(bowTie); + + expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(bowTie, referenceSquare), 0); + + // A reference in the empty wedge below the crossing point is not intersected. + const inWedge = { type: "Polygon" as const, coordinates: [rectangle(2.98, 48.41, 3.02, 48.43)] }; + expect(deriveExtras(bowTie, inWedge).intersection_area).toEqual(0); + + // Nor is a bow tie lying inside the reference reduced to its signed area. + const inside = { type: "Polygon" as const, coordinates: [rectangle(2, 48, 4, 49)] }; + expect(deriveExtras(bowTie, inside).intersection_area).toBeCloseTo(polyclipIntersectionArea(bowTie, inside), 0); + }); + it("should return null for a non-areal geometry", () => { const feature = deriveExtras({ type: "LineString", coordinates: [[2.5, 48.5], [3.5, 48.5]] }); From 8c275db125c0334fe2474ba61ff52ecaae86a495 Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:47:11 +0200 Subject: [PATCH 5/6] test(intersection_area): add a benchmark on synthetic geometries npm run bench runs test/wfs/intersectionArea.bench.ts (vitest bench) through transformFeatureCollectionResponse, as for one request: small polygons inside and across a detailed boundary, large polygons with holes across it, and sectors sharing it position for position. --- package.json | 1 + test/wfs/intersectionArea.bench.ts | 78 ++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+) create mode 100644 test/wfs/intersectionArea.bench.ts diff --git a/package.json b/package.json index f0450efb..41fa89a4 100644 --- a/package.json +++ b/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", diff --git a/test/wfs/intersectionArea.bench.ts b/test/wfs/intersectionArea.bench.ts new file mode 100644 index 00000000..18843f4e --- /dev/null +++ b/test/wfs/intersectionArea.bench.ts @@ -0,0 +1,78 @@ +/** + * Performance of `intersection_area`, through `transformFeatureCollectionResponse` as for + * one request, on synthetic geometries reproducing the costly cases met on real data. + * + * npm run bench + * + * To compare two revisions: `npm run bench -- --outputJson before.json` on the first one, + * then `npm run bench -- --compare before.json` on the second one. + */ +import { bench, describe } from "vitest"; +import type { Polygon, Position } from "geojson"; + +import { transformFeatureCollectionResponse } from "../../src/wfs/response.js"; + +const GOLDEN_ANGLE = Math.PI * (3 - Math.sqrt(5)); + +/** Closed ring of `vertices` positions around a center, its radius waving `waves` times. */ +function ring(lon: number, lat: number, radius: number, vertices: number, waves = 0) : Position[] { + const positions = Array.from({ length: vertices }, (_, index) => { + const angle = (Math.PI * 2 * index) / vertices; + const r = radius * (1 + 0.02 * Math.sin(waves * angle)); + return [lon + r * Math.cos(angle), lat + r * Math.sin(angle)]; + }); + return [...positions, positions[0]]; +} + +function polygon(...rings: Position[][]) : Polygon { + return { type: "Polygon", coordinates: rings }; +} + +// A region with a detailed boundary, like a département (about 60,000 positions in BD TOPO). +const [LON, LAT, RADIUS] = [2.35, 48.85, 0.3]; +const regionRing = ring(LON, LAT, RADIUS, 20000, 200); +const region = polygon(regionRing); + +/** Position on the region boundary at the given angle. */ +function boundary(angle: number) : [number, number] { + const r = RADIUS * (1 + 0.02 * Math.sin(200 * angle)); + return [LON + r * Math.cos(angle), LAT + r * Math.sin(angle)]; +} + +const scenarios: Record = { + // Vegetation zones inside a département: spread over a disc, well inside the boundary. + "5000 small polygons inside a detailed region": Array.from({ length: 5000 }, (_, index) => { + const r = 0.25 * Math.sqrt((index + 0.5) / 5000), angle = index * GOLDEN_ANGLE; + return polygon(ring(LON + r * Math.cos(angle), LAT + r * Math.sin(angle), 0.001, 12)); + }), + // Every polygon needs an actual intersection. + "5000 small polygons across a detailed boundary": Array.from({ length: 5000 }, (_, index) => + polygon(ring(...boundary((Math.PI * 2 * index) / 5000), 0.002, 12))), + // Large polygons with holes across the boundary: both sides are detailed. + "20 large polygons with holes across a detailed boundary": Array.from({ length: 20 }, (_, index) => { + const [lon, lat] = boundary((Math.PI * 2 * index) / 20); + const holes = [-0.03, 0.03].map((offset) => ring(lon + offset, lat + offset, 0.01, 256).reverse()); + return polygon(ring(lon, lat, 0.1, 20000, 100), ...holes); + }), + // Départements of a région: pie slices sharing the region boundary, position for position. + "20 sectors tiling the region, sharing its boundary": Array.from({ length: 20 }, (_, index) => { + const arc = regionRing.slice(index * 1000, (index + 1) * 1000 + 1); + return polygon([[LON, LAT], ...arc, [LON, LAT]]); + }), +}; + +for (const [name, geometries] of Object.entries(scenarios)) { + const featureCollection = { + type: "FeatureCollection", + features: geometries.map((geometry, index) => ({ id: `feature.${index}`, geometry, properties: {} })), + }; + describe(name, () => { + bench("intersection_area", () => { + transformFeatureCollectionResponse(featureCollection, { + typename: "BENCH:feature", + spatial_extras: ["intersection_area"], + intersects_feature_filter: { typename: "BENCH:region", feature_id: "region.1" }, + }, region); + }); + }); +} From 0655b8b564a10bf2840817abc83c2907ae3e1417 Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Wed, 30 Sep 2026 16:38:42 +0200 Subject: [PATCH 6/6] perf: return null instead of forwarding to turf on invalid geometries --- package-lock.json | 8 +- package.json | 2 +- src/helpers/intersectionArea.ts | 186 ++++++++++---------------- test/helpers/intersectionArea.test.ts | 34 ++++- test/wfs/response.test.ts | 21 +-- 5 files changed, 121 insertions(+), 130 deletions(-) diff --git a/package-lock.json b/package-lock.json index a4c25ccf..cf637ac6 100644 --- a/package-lock.json +++ b/package-lock.json @@ -18,7 +18,6 @@ "@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", @@ -40,6 +39,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", @@ -1806,6 +1806,7 @@ "version": "7.4.0", "resolved": "https://registry.npmjs.org/@turf/intersect/-/intersect-7.4.0.tgz", "integrity": "sha512-n/wbZfDPoM1JKwrql8mHBoS9cdRWUC7/VE9IBZOttVPv8Guo0B6+Mzc2Z2e1PdIkSP/1R2yH6lxPpWMuuBfVEA==", + "dev": true, "license": "MIT", "dependencies": { "@turf/helpers": "7.4.0", @@ -1822,6 +1823,7 @@ "version": "7.4.0", "resolved": "https://registry.npmjs.org/@turf/helpers/-/helpers-7.4.0.tgz", "integrity": "sha512-7PAwLZqOdRzTI5g9bHvUwlloAXPDH/mlajtryk0tw4ZwGMtmXsAyF4QundsAMfy4u48Fyd3AUqMXY0SMxqJGWg==", + "dev": true, "license": "MIT", "dependencies": { "@types/geojson": "^7946.0.10", @@ -1835,6 +1837,7 @@ "version": "7.4.0", "resolved": "https://registry.npmjs.org/@turf/meta/-/meta-7.4.0.tgz", "integrity": "sha512-3cLUvlEyDuSnMSzrjhaLAEiYR8xhbfWyTVQlrlEw40xL81d4KF4PqUWbjTXKpXZStdYbet2GCurl79KMypNw6g==", + "dev": true, "license": "MIT", "dependencies": { "@turf/helpers": "7.4.0", @@ -2772,6 +2775,7 @@ "version": "9.3.1", "resolved": "https://registry.npmjs.org/bignumber.js/-/bignumber.js-9.3.1.tgz", "integrity": "sha512-Ko0uX15oIUS7wJ3Rb30Fs6SkVbLmPBAKdlm7q9+ak9bbIeFf0MwuBsQV6z7+X768/cHsfg+WlysDWJcmthjsjQ==", + "dev": true, "license": "MIT", "engines": { "node": "*" @@ -5546,6 +5550,7 @@ "version": "0.16.8", "resolved": "https://registry.npmjs.org/polyclip-ts/-/polyclip-ts-0.16.8.tgz", "integrity": "sha512-JPtKbDRuPEuAjuTdhR62Gph7Is2BS1Szx69CFOO3g71lpJDFo78k4tFyi+qFOMVPePEzdSKkpGU3NBXPHHjvKQ==", + "dev": true, "license": "MIT", "dependencies": { "bignumber.js": "^9.1.0", @@ -6077,6 +6082,7 @@ "version": "1.0.2", "resolved": "https://registry.npmjs.org/splaytree-ts/-/splaytree-ts-1.0.2.tgz", "integrity": "sha512-0kGecIZNIReCSiznK3uheYB8sbstLjCZLiwcQwbmLhgHJj2gz6OnSPkVzJQCMnmEz1BQ4gPK59ylhBoEWOhGNA==", + "dev": true, "license": "BDS-3-Clause" }, "node_modules/split2": { diff --git a/package.json b/package.json index 41fa89a4..9e588792 100644 --- a/package.json +++ b/package.json @@ -64,7 +64,6 @@ "@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", @@ -82,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", diff --git a/src/helpers/intersectionArea.ts b/src/helpers/intersectionArea.ts index 2fb270cf..63deac42 100644 --- a/src/helpers/intersectionArea.ts +++ b/src/helpers/intersectionArea.ts @@ -1,11 +1,7 @@ import earcut, { deviation } from "earcut"; import { bbox } from "@turf/bbox"; -import { bboxClip } from "@turf/bbox-clip"; -import { intersect } from "@turf/intersect"; -import { feature, featureCollection } from "@turf/helpers"; import type { BBox, MultiPolygon, Polygon, Position } from "geojson"; import area, { EARTH_RADIUS } from "./area.js"; -import { dropEmptyRings } from "./geojson.js"; /** * Area of the part of areal features lying inside an areal filter. @@ -25,8 +21,10 @@ import { dropEmptyRings } from "./geojson.js"; * the product of both multiplicities: the polygons of a MultiPolygon are summed, overlaps * included, as `area` does. * - * Self-intersecting rings (a bow tie, a clipped piece with the wrong orientation, a proper - * crossing where earcut fails) and oversized leaves are left to `@turf/intersect`. + * A ring detected as self-intersecting (no area but some extent, like a bow tie; a clipped + * piece with the wrong orientation; a proper crossing where earcut fails) throws an + * `InvalidGeometryError`. Detection is not exhaustive: a self-intersecting ring whose pieces + * all keep its orientation counts with its winding numbers, as in `area`. * * Areas are exact on the sphere for edges straight in lon/lat (same integral and radius as * `area`), evaluated relative to a local origin to avoid cancellation on small rings. @@ -35,13 +33,15 @@ import { dropEmptyRings } from "./geojson.js"; /** Area (m²) of the part of an areal feature lying inside the prepared filter. */ export type IntersectionArea = (geo: Polygon | MultiPolygon) => number; -/** A box becomes a leaf when small × big <= max(LEAF_PRODUCT, LEAF_FACTOR × (small + big)) positions. */ +/** + * A box becomes a leaf when small × big <= min(MAX_LEAF_OPERATIONS, max(LEAF_PRODUCT, LEAF_FACTOR × (small + big))) + * positions, a leaf costing about small × big triangle × position operations. + */ const LEAF_PRODUCT = 1024; const LEAF_FACTOR = 16; -/** Bound on the recursion, in case positions pile up at the same place. */ -const MAX_DEPTH = 60; -/** Triangle × position operations of a leaf above which `@turf/intersect` takes over. */ const MAX_LEAF_OPERATIONS = 1e6; +/** Bound on the recursion, in case positions pile up at the same place: the leaf is then computed whatever its size. */ +const MAX_DEPTH = 60; /** Relative snap to the feature area, and to 0 (tighter, so that real slivers survive). */ const SNAP_FULL = 1e-9; const SNAP_ZERO = 1e-12; @@ -52,8 +52,13 @@ const DEGREES_TO_RADIANS = Math.PI / 180; const HALF_DEGREES_TO_RADIANS = DEGREES_TO_RADIANS / 2; const R2 = EARTH_RADIUS ** 2; -/** Thrown to hand a feature over to `@turf/intersect`. */ -const FALLBACK = Symbol("fallback"); +/** Thrown when a ring of the feature or of the filter is detected as self-intersecting. */ +export class InvalidGeometryError extends Error { + constructor(message = "Un anneau auto-intersectant a été détecté dans l'objet ou dans le filtre spatial.") { + super(message); + this.name = "InvalidGeometryError"; + } +} /** * A ring in flat local coordinates `[x0, y0, x1, y1, ...]`, without its closing position, @@ -110,30 +115,33 @@ function makeRing(src: Float64Array, count2: number, mu: number, o: number) : Ri // rounding, reveals a self-intersecting ring. const extent = Math.max(x1 - x0, y1 - y0); const tolerance = 2e-6 * (x1 - x0) * (y1 - y0) + 1e-12 * n * extent * extent; - if (area2 * o < -tolerance) throw FALLBACK; + if (area2 * o < -tolerance) throw new InvalidGeometryError(); return { c, n, x0, y0, x1, y1, mu, o }; } -/** Local copy of a GeoJSON ring (x = lon - lon0, y = lat - lat0) without its closing position, or null below 3 positions. */ +/** + * Local copy of a GeoJSON ring (x = lon - lon0, y = lat - lat0), or null below 3 positions. + * Repeated consecutive positions, the closing one included, are dropped: they add no area, + * only work, and would pile up in the smallest boxes. + */ function localRing(positions: Position[], lon0: number, lat0: number) : Ring | null { - let n = positions.length; - if (n > 1) { - const first = positions[0], last = positions[n - 1]; - if (first[0] === last[0] && first[1] === last[1]) n--; - } - if (n < 3) return null; - const c = new Float64Array(2 * n); + const c = new Float64Array(2 * positions.length); + let n = 0; let x0 = Infinity, y0 = Infinity, x1 = -Infinity, y1 = -Infinity; - for (let i = 0; i < n; i++) { - const x = positions[i][0] - lon0, y = positions[i][1] - lat0; - c[2 * i] = x; - c[2 * i + 1] = y; + for (const position of positions) { + const x = position[0] - lon0, y = position[1] - lat0; + if (n > 0 && x === c[2 * n - 2] && y === c[2 * n - 1]) continue; + c[2 * n] = x; + c[2 * n + 1] = y; + n++; if (x < x0) x0 = x; if (x > x1) x1 = x; if (y < y0) y0 = y; if (y > y1) y1 = y; } - return { c, n, x0, y0, x1, y1, mu: 0, o: 0 }; + while (n > 1 && c[2 * n - 2] === c[0] && c[2 * n - 1] === c[1]) n--; + if (n < 3) return null; + return { c: c.subarray(0, 2 * n), n, x0, y0, x1, y1, mu: 0, o: 0 }; } /** Twice the planar signed area (positive counterclockwise), as a fan from the first position. */ @@ -383,7 +391,7 @@ function boxIntegral(F: Ring[], G: Ring[], x0: number, y0: number, x1: number, y function solve(F: Ring[], G: Ring[], x0: number, y0: number, x1: number, y1: number, depth: number, ctx: AreaContext) : number { const nF = countPositions(F), nG = countPositions(G); const small = Math.min(nF, nG), big = Math.max(nF, nG); - if (depth >= MAX_DEPTH || small * big <= Math.max(LEAF_PRODUCT, LEAF_FACTOR * (small + big))) { + if (depth >= MAX_DEPTH || small * big <= Math.min(MAX_LEAF_OPERATIONS, Math.max(LEAF_PRODUCT, LEAF_FACTOR * (small + big)))) { return leaf(F, G, nF, nG, ctx); } const FL: Ring[] = [], FH: Ring[] = [], GL: Ring[] = [], GH: Ring[] = []; @@ -404,9 +412,7 @@ function solve(F: Ring[], G: Ring[], x0: number, y0: number, x1: number, y1: num // --- Leaves: triangulate the smaller side, clip the other side's rings by each triangle --- function leaf(F: Ring[], G: Ring[], nF: number, nG: number, ctx: AreaContext) : number { - const [A, B, nA, nB]: [Ring[], Ring[], number, number] = nF <= nG ? [F, G, nF, nG] : [G, F, nG, nF]; - // About one triangle per position of A, each clipping every position of B. - if (nA * nB > MAX_LEAF_OPERATIONS) throw FALLBACK; + const [A, B] = nF <= nG ? [F, G] : [G, F]; let total = 0; for (const a of A) total += a.mu * ringIntegral(a, B, ctx); return total; @@ -429,8 +435,8 @@ function ringIntegral(a: Ring, B: Ring[], ctx: AreaContext) : number { return area2 > 0 ? total : -total; } // earcut fails on the zero-width bridges left by Sutherland-Hodgman, which a signed fan - // handles exactly, but also on self-intersecting rings, left to `@turf/intersect`. - if (hasProperCrossing(c, n)) throw FALLBACK; + // handles exactly, but also on self-intersecting rings, rejected. + if (hasProperCrossing(c, n)) throw new InvalidGeometryError(); // Signed fan from the first position: winding(a) = sum of the signed windings of its triangles. let total = 0; const ax = c[0], ay = c[1]; @@ -494,7 +500,7 @@ function triangleIntegral(B: Ring[], ax: number, ay: number, bx: number, by: num for (const r of B) { if (r.x1 < tx0 || r.x0 > tx1 || r.y1 < ty0 || r.y0 > ty1) continue; const clipped = clippedArea(r, ax, ay, bx, by, cx, cy, ctx); - if (clipped * r.o < -(relative + rounding * r.n)) throw FALLBACK; + if (clipped * r.o < -(relative + rounding * r.n)) throw new InvalidGeometryError(); total += r.mu * clipped; } return total; @@ -557,8 +563,6 @@ type PreparedFilter = { /** Local origin: the centre of the filter bbox. */ ox: number; oy: number; - /** True when a filter ring has no area but some extent (self-intersecting, like a bow tie). */ - invalid: boolean; /** Rings of the tile containing an absolute bbox, in local coordinates. */ near: (minX: number, minY: number, maxX: number, maxY: number) => Ring[]; }; @@ -574,14 +578,14 @@ function prepareFilter(filter: Polygon | MultiPolygon) : PreparedFilter { const ox = (west + east) / 2, oy = (south + north) / 2; const size = Math.max(east - west, north - south); const root: Ring[] = []; - let invalid = false; for (const polygon of polygonsOf(filter)) { for (let index = 0; index < polygon.length; index++) { const ring = localRing(polygon[index], ox, oy); if (!ring) continue; const area2 = planarArea2(ring.c, ring.n); if (Math.abs(area2) <= 1e-12 * 2 * (ring.x1 - ring.x0) * (ring.y1 - ring.y0)) { - if (ring.x1 > ring.x0 && ring.y1 > ring.y0) invalid = true; + // No area but some extent: self-intersecting (a bow tie) or zero-width. Without extent, it covers nothing. + if (ring.x1 > ring.x0 && ring.y1 > ring.y0) throw new InvalidGeometryError("Le filtre spatial contient un anneau auto-intersectant."); continue; } ring.o = area2 > 0 ? 1 : -1; @@ -604,7 +608,7 @@ function prepareFilter(filter: Polygon | MultiPolygon) : PreparedFilter { return rings; } return { - ox, oy, invalid, + ox, oy, near(minX, minY, maxX, maxY) { // Deepest tile whose cell is at least twice as wide as the bbox. const width = Math.max(maxX - minX, maxY - minY); @@ -620,25 +624,18 @@ function prepareFilter(filter: Polygon | MultiPolygon) : PreparedFilter { /** * Prepare the filter once (per request), and return the function computing the area (m²) - * of the part of each feature lying inside it. + * of the part of each feature lying inside it. Both throw an `InvalidGeometryError` on a + * ring detected as self-intersecting. */ export function makeIntersectionArea(filter: Polygon | MultiPolygon) : IntersectionArea { const prepared = prepareFilter(filter); return (geo) => { const featureBbox = bbox(geo); - if (prepared.invalid) return filterFallback(filter, geo, featureBbox); // The filter near the feature: its tile, then exactly the feature bbox. const [f0, f1, f2, f3] = featureBbox; const { ox, oy } = prepared; const fx0 = f0 - ox, fy0 = f1 - oy, fx1 = f2 - ox, fy1 = f3 - oy; - let near: Ring[]; - try { - near = clipRings(prepared.near(f0, f1, f2, f3), fx0, fy0, fx1, fy1); - } catch (error) { - // A tile piece with the wrong orientation: self-intersecting filter ring. - if (error !== FALLBACK) throw error; - return filterFallback(filter, geo, featureBbox); - } + const near = clipRings(prepared.near(f0, f1, f2, f3), fx0, fy0, fx1, fy1); if (near.length == 0) return 0; let W0 = Infinity, W1 = Infinity, W2 = -Infinity, W3 = -Infinity; for (const r of near) { @@ -654,37 +651,14 @@ export function makeIntersectionArea(filter: Polygon | MultiPolygon) : Intersect const result = Number.isNaN(constant) ? null : insideShortcut(constant, geo, featureBbox); if (result !== null) return result; } - const fallbackInputs = () : [Polygon | MultiPolygon, BBox] | null => { - const clippedFilter = dropEmptyRings(bboxClip(filter, featureBbox).geometry); - return clippedFilter ? [clippedFilter, bbox(clippedFilter)] : null; - }; - return featureIntegral(geo, near, W0, W1, W2, W3, ox, oy, fallbackInputs); + return featureIntegral(geo, near, W0, W1, W2, W3, ox, oy); }; } -/** `@turf/intersect` on the feature and the whole filter clipped to its bbox. */ -function filterFallback(filter: Polygon | MultiPolygon, geo: Polygon | MultiPolygon, featureBbox: BBox) : number { - const clippedFilter = dropEmptyRings(bboxClip(filter, featureBbox).geometry); - if (!clippedFilter) return 0; - return snap(polyclipFallback(geo, clippedFilter, bbox(clippedFilter)), geo, area(geo), false); -} - -/** `@turf/intersect` on the pre-clipped pair, one feature polygon at a time, so that polygons are summed. */ -function polyclipFallback(geo: Polygon | MultiPolygon, clippedFilter: Polygon | MultiPolygon, workingBox: BBox) : number { - let total = 0; - for (const polygon of polygonsOf(geo)) { - const clipped = dropEmptyRings(bboxClip({ type: "Polygon", coordinates: polygon }, workingBox).geometry); - if (!clipped) continue; - const inter = intersect(featureCollection([feature(clipped), feature(clippedFilter)])); - total += inter == null ? 0 : area(inter.geometry); - } - return total; -} - /** * Result when the filter has a constant multiplicity over the whole feature bbox, or null to * take the general path: a feature without area but with some extent may be a bow tie, - * whose `area` is not the area it covers. + * which the general path rejects. */ function insideShortcut(constant: number, geo: Polygon | MultiPolygon, [f0, f1, f2, f3]: BBox) : number | null { if (constant <= 0) return 0; @@ -703,7 +677,6 @@ function featureIntegral( G: Ring[], bx0: number, by0: number, bx1: number, by1: number, lon0: number, lat0: number, - fallbackInputs: () => [Polygon | MultiPolygon, BBox] | null, ) : number { const phi0 = lat0 * DEGREES_TO_RADIANS; const ctx: AreaContext = { phi0, s0: Math.sin(phi0), c0: Math.cos(phi0) }; @@ -712,7 +685,6 @@ function featureIntegral( const F: Ring[] = []; const whole: Ring[] = []; // every ring of the feature, unclipped, for the snapping decision let clipped = false; - let invalid = false; let planarFeature2 = 0; // twice the planar area of the feature, for the snapping estimate for (const polygon of polygonsOf(geo)) { for (let index = 0; index < polygon.length; index++) { @@ -721,9 +693,9 @@ function featureIntegral( const area2 = planarArea2(ring.c, ring.n); const bboxArea2 = 2 * (ring.x1 - ring.x0) * (ring.y1 - ring.y0); if (Math.abs(area2) <= 1e-12 * bboxArea2) { - // No area but some extent: a self-intersecting ring whose lobes cancel (bow tie) - // is left to `@turf/intersect`; a zero-width ring is only slower there. - if (bboxArea2 > 0) invalid = true; + // No area but some extent: self-intersecting (a bow tie whose lobes cancel) or + // zero-width. Without extent, it covers nothing. + if (bboxArea2 > 0) throw new InvalidGeometryError("L'objet contient un anneau auto-intersectant."); continue; } const role = index == 0 ? 1 : -1; @@ -736,57 +708,41 @@ function featureIntegral( continue; } clipped = true; - if (ring.x1 < bx0 || ring.x0 > bx1 || ring.y1 < by0 || ring.y0 > by1 || invalid) continue; - try { - const piece = clipRingToBox(ring, bx0, by0, bx1, by1); - if (piece) F.push(piece); - } catch (error) { - if (error !== FALLBACK) throw error; - invalid = true; // a piece with the wrong orientation: self-intersecting ring - } + if (ring.x1 < bx0 || ring.x0 > bx1 || ring.y1 < by0 || ring.y0 > by1) continue; + const piece = clipRingToBox(ring, bx0, by0, bx1, by1); + if (piece) F.push(piece); } } const estimate = planarFeature2 / 2 * DEGREES_TO_RADIANS * DEGREES_TO_RADIANS * ctx.c0 * R2; - if (!invalid) { - try { - const F2: Ring[] = [], G2: Ring[] = []; - const cF = extractConstant(F, bx0, by0, bx1, by1, F2); - const cG = extractConstant(G, bx0, by0, bx1, by1, G2); - if (G2.length == 0 && !clipped) { - // The filter covers the whole working box, which contains the feature, or none of it. - return cG > 0 ? area(geo) : 0; - } - let total = 0; - if (cF !== 0 && cG !== 0) total += cF * cG * boxArea(bx0, by0, bx1, by1, ctx); - if (cF !== 0 && G2.length) total += cF * sumArea(G2, ctx); - if (cG !== 0 && F2.length) total += cG * sumArea(F2, ctx); - if (F2.length && G2.length) total += solve(F2, G2, bx0, by0, bx1, by1, 0, ctx); - // Capped at the feature area. - return snap(total * R2, geo, estimate, true, () => sumArea(whole, ctx) * R2); - } catch (error) { - if (error !== FALLBACK) throw error; - } + const F2: Ring[] = [], G2: Ring[] = []; + const cF = extractConstant(F, bx0, by0, bx1, by1, F2); + const cG = extractConstant(G, bx0, by0, bx1, by1, G2); + if (G2.length == 0 && !clipped) { + // The filter covers the whole working box, which contains the feature, or none of it. + return cG > 0 ? area(geo) : 0; } - // Not capped: for a self-intersecting feature, `area` is not the area it covers. - const inputs = fallbackInputs(); - if (!inputs) return 0; - return snap(polyclipFallback(geo, inputs[0], inputs[1]), geo, estimate, false); + let total = 0; + if (cF !== 0 && cG !== 0) total += cF * cG * boxArea(bx0, by0, bx1, by1, ctx); + if (cF !== 0 && G2.length) total += cF * sumArea(G2, ctx); + if (cG !== 0 && F2.length) total += cG * sumArea(F2, ctx); + if (F2.length && G2.length) total += solve(F2, G2, bx0, by0, bx1, by1, 0, ctx); + return snap(total * R2, geo, estimate, () => sumArea(whole, ctx) * R2); } /** - * Snap to the feature area (within SNAP_FULL, or above it when `cap`) or to 0 (within + * Snap to the feature area (within SNAP_FULL, or above it: capped) or to 0 (within * SNAP_ZERO). The exact feature area is only computed when the planar estimate says that * the result may be close to either. */ -function snap(total: number, geo: Polygon | MultiPolygon, estimate: number, cap: boolean, localFeatureArea?: () => number) : number { +function snap(total: number, geo: Polygon | MultiPolygon, estimate: number, localFeatureArea: () => number) : number { if (Number.isNaN(total)) throw new Error("intersection area is NaN"); if (total <= 0) return 0; if (total <= 100 * SNAP_ZERO * estimate || total >= 0.5 * estimate) { - // Decided on the same (local) formula as `total` when available; the snapped value is - // `area`, so that a feature inside the filter gets exactly its `area`. - const featureArea = localFeatureArea ? localFeatureArea() : area(geo); - if (total >= featureArea * (1 - SNAP_FULL) && (cap || total <= featureArea * (1 + SNAP_FULL))) return area(geo); + // Decided on the same (local) formula as `total`; the snapped value is `area`, so that + // a feature inside the filter gets exactly its `area`. + const featureArea = localFeatureArea(); + if (total >= featureArea * (1 - SNAP_FULL)) return area(geo); if (total <= featureArea * SNAP_ZERO) return 0; } return total; diff --git a/test/helpers/intersectionArea.test.ts b/test/helpers/intersectionArea.test.ts index 5579a202..966eff06 100644 --- a/test/helpers/intersectionArea.test.ts +++ b/test/helpers/intersectionArea.test.ts @@ -4,7 +4,7 @@ import { feature, featureCollection } from "@turf/helpers"; import type { MultiPolygon, Polygon, Position } from "geojson"; import area from "../../src/helpers/area.js"; -import { makeIntersectionArea } from "../../src/helpers/intersectionArea.js"; +import { InvalidGeometryError, makeIntersectionArea } from "../../src/helpers/intersectionArea.js"; function rectangle(west: number, south: number, east: number, north: number) : Position[] { return [[west, south], [east, south], [east, north], [west, north], [west, south]]; @@ -57,13 +57,39 @@ describe("helpers/intersectionArea", () => { } }); - it("should fall back on general polygon clipping for a self-intersecting filter", () => { + it("should ignore repeated positions", () => { + // Both boundaries pass through the tip of a notch, repeated on both sides: the + // repetitions all end in the deepest boxes around it. + const tip: Position = [2.95, 48.5]; + const withTip = (times: number) => Array.from({ length: times }, () => tip); + const notchedFilter = (times: number) : Polygon => ({ type: "Polygon", coordinates: [[[2, 48], [3, 48], ...withTip(times), [3, 49], [2, 49], [2, 48]]] }); + const notchedGeo = (times: number) : Polygon => ({ type: "Polygon", coordinates: [[[2.9, 48.45], [3.05, 48.45], ...withTip(times), [3.05, 48.55], [2.9, 48.55], [2.9, 48.45]]] }); + + expect(makeIntersectionArea(notchedFilter(100000))(notchedGeo(100000))).toEqual(makeIntersectionArea(notchedFilter(1))(notchedGeo(1))); + }); + + it("should reject a self-intersecting filter", () => { const bowTie: Polygon = { type: "Polygon", coordinates: [[[2, 48], [3, 49], [3, 48], [2, 49], [2, 48]]], }; - const geo: Polygon = { type: "Polygon", coordinates: [rectangle(2.1, 48.1, 2.9, 48.4)] }; - expect(makeIntersectionArea(bowTie)(geo)).toBeCloseTo(polyclipArea(geo, bowTie), 0); + expect(() => makeIntersectionArea(bowTie)).toThrow(InvalidGeometryError); + }); + + it("should reject a self-intersecting geometry whose part inside the filter has the wrong orientation", () => { + // Figure eight: its larger lobe, counterclockwise, lies mostly east of the filter, its + // smaller lobe, clockwise, inside it, so that the piece west of the eastern edge is clockwise. + const larger = Array.from({ length: 40 }, (_, index) => { + const angle = Math.PI + (Math.PI * 2 * index) / 40; + return [3.08 + 0.1 * Math.cos(angle), 48.5 + 0.1 * Math.sin(angle)]; + }); + const smaller = Array.from({ length: 40 }, (_, index) => { + const angle = (Math.PI * 2 * index) / 40; + return [2.95 + 0.03 * Math.cos(angle), 48.5 - 0.03 * Math.sin(angle)]; + }); + const figureEight: Polygon = { type: "Polygon", coordinates: [[...larger, ...smaller, larger[0]]] }; + + expect(() => makeIntersectionArea(filter)(figureEight)).toThrow(InvalidGeometryError); }); }); diff --git a/test/wfs/response.test.ts b/test/wfs/response.test.ts index 113fcd1b..ab0599e4 100644 --- a/test/wfs/response.test.ts +++ b/test/wfs/response.test.ts @@ -670,23 +670,26 @@ describe("wfs_engine/response", () => { expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(withHole, referenceSquare), 0); }); - it("should fall back on general polygon clipping for a self-intersecting geometry", () => { + it("should return null for a self-intersecting geometry", () => { // Bow tie, whose lobes cancel in its signed area. const bowTie: Polygon = { type: "Polygon", coordinates: [[[2.8, 48.4], [3.2, 48.6], [3.2, 48.4], [2.8, 48.6], [2.8, 48.4]]], }; - const feature = deriveExtras(bowTie); - expect(feature.intersection_area).toBeCloseTo(polyclipIntersectionArea(bowTie, referenceSquare), 0); - - // A reference in the empty wedge below the crossing point is not intersected. + // Across the reference boundary, around its crossing point, and inside the reference. + expect(deriveExtras(bowTie).intersection_area).toBeNull(); const inWedge = { type: "Polygon" as const, coordinates: [rectangle(2.98, 48.41, 3.02, 48.43)] }; - expect(deriveExtras(bowTie, inWedge).intersection_area).toEqual(0); - - // Nor is a bow tie lying inside the reference reduced to its signed area. + expect(deriveExtras(bowTie, inWedge).intersection_area).toBeNull(); const inside = { type: "Polygon" as const, coordinates: [rectangle(2, 48, 4, 49)] }; - expect(deriveExtras(bowTie, inside).intersection_area).toBeCloseTo(polyclipIntersectionArea(bowTie, inside), 0); + expect(deriveExtras(bowTie, inside).intersection_area).toBeNull(); + }); + + it("should return null for a self-intersecting reference", () => { + const bowTie = { type: "Polygon" as const, coordinates: [[[2, 48], [3, 49], [3, 48], [2, 49], [2, 48]]] }; + const feature = deriveExtras({ type: "Polygon", coordinates: [rectangle(2.1, 48.1, 2.9, 48.4)] }, bowTie); + + expect(feature.intersection_area).toBeNull(); }); it("should return null for a non-areal geometry", () => {