/** * Spatial analysis library for WGS84 coordinates using spherical Earth model * Supports: Point, LineString, Polygon, MultiLineString * @module Geometry */ // Earth radius in meters (simplified sphere model) const R = 6371000; const TOLERANCE = 1e-10; /** * Geographic bounds representation */ export class Bounds { /** * @param {number} minX - Minimum longitude * @param {number} minY - Minimum latitude * @param {number} maxX - Maximum longitude * @param {number} maxY - Maximum latitude */ constructor(minX, minY, maxX, maxY) { this.minX = minX; this.minY = minY; this.maxX = maxX; this.maxY = maxY; } /** * Create a union of multiple Bounds objects to cover the combined area * @param {Bounds[]} boundsArray - Array of Bounds instances * @returns {Bounds} A new Bounds instance covering all input bounds */ static union(boundsArray) { const minX = Math.min(...boundsArray.map((b) => b.minX)); const minY = Math.min(...boundsArray.map((b) => b.minY)); const maxX = Math.max(...boundsArray.map((b) => b.maxX)); const maxY = Math.max(...boundsArray.map((b) => b.maxY)); return new Bounds(minX, minY, maxX, maxY); } /** @returns {number} Width in degrees */ get width() { return this.maxX - this.minX; } /** @returns {number} Height in degrees */ get height() { return this.maxY - this.minY; } /** @returns {Coordinate} Center coordinate */ get center() { return [this.minX + this.width / 2, this.minY + this.height / 2]; } /** * Check if bounds contain coordinate * @param {Coordinate} coordinate * @returns {boolean} */ contains(coordinate) { const [lng, lat] = coordinate; return lng >= this.minX && lng <= this.maxX && lat >= this.minY && lat <= this.maxY; } } /** * @typedef {[number, number]} Coordinate - [longitude, latitude] in decimal degrees */ /** * Base geometry class * @abstract */ export class Geometry { /** @returns {string} Geometry type */ get type() { return this.constructor.name; } /** * Convert to GeoJSON geometry object * @returns {Object} GeoJSON geometry * @abstract */ toGeoJSON() { throw new Error("Not implemented"); } /** * Calculate geometry bounds * @returns {Bounds} Geometry bounds * @abstract */ bounds() { throw new Error("Not implemented"); } /** * Calculate distance to another geometry * @param {Geometry} other - Target geometry * @returns {number} Distance in meters */ distance(other) { return Geometry.distance(this, other); } /** * Calculate geometry length * @returns {number} Length in meters */ length() { return Geometry.length(this); } /** * Check if geometry intersects another * @param {Geometry} other - Target geometry * @returns {boolean} */ intersects(other) { return Geometry.intersects(this, other); } /** * Check if geometry is completely within another * @param {Geometry} other - Container geometry * @returns {boolean} */ within(other) { return Geometry.within(this, other); } /** * Deserialize geometry from GeoJSON * @param {{type: string, coordinates?: any, [key: string]: any}} geojson * @returns {Geometry|null} */ static fromGeoJSON(geojson) { if (!geojson?.type) return null; switch (geojson.type) { case "Point": return Point.fromGeoJSON(geojson); case "LineString": return LineString.fromGeoJSON(geojson); case "Polygon": return Polygon.fromGeoJSON(geojson); case "MultiLineString": return MultiLineString.fromGeoJSON(geojson); default: throw new Error(`Invalid GeoJSON object: missing required 'type' property`); } } /** * Compute geometry bounds * @param {Geometry} geometry - Input geometry * @returns {Bounds} Geometry bounds */ static bounds(geometry) { if (geometry instanceof Point) { return new Bounds(geometry.lng, geometry.lat, geometry.lng, geometry.lat); } if (geometry instanceof LineString) { return Geometry.calculateBoundsFromCoords(geometry.coordinates); } if (geometry instanceof Polygon) { return Geometry.calculateBoundsFromCoords(geometry.rings[0]); } throw new Error(`Unsupported geometry type: ${geometry.type}`); } /** * Calculate bounds from coordinate array * @param {Coordinate[]} coords - Coordinate array * @returns {Bounds} Calculated bounds */ static calculateBoundsFromCoords(coords) { const lngs = coords.map(([lng]) => lng); const lats = coords.map(([, lat]) => lat); return new Bounds(Math.min(...lngs), Math.min(...lats), Math.max(...lngs), Math.max(...lats)); } /** * Calculate distance between geometries * @param {Geometry} a - First geometry * @param {Geometry} b - Second geometry * @returns {number} Minimum distance in meters */ static distance(a, b) { const pointsA = Geometry.#geometryToPoints(a); const pointsB = Geometry.#geometryToPoints(b); let minDistance = Infinity; for (const pointA of pointsA) { for (const pointB of pointsB) { const dist = Point.distance(pointA, pointB); if (dist < minDistance) minDistance = dist; } } return minDistance; } /** * Convert geometry to representative points * @param {Geometry} geometry - Input geometry * @returns {Point[]} Array of points */ static #geometryToPoints(geometry) { if (geometry instanceof Point) { return [geometry]; } if (geometry instanceof LineString) { return geometry.coordinates.map(([lng, lat]) => new Point(lng, lat)); } if (geometry instanceof Polygon) { return geometry.rings.flatMap((ring) => ring.map(([lng, lat]) => new Point(lng, lat))); } throw new Error(`Unsupported geometry type: ${geometry.type}`); } /** * Calculate geometry length * @param {Geometry} geometry - Input geometry * @returns {number} Length in meters */ static length(geometry) { if (geometry instanceof Point) return 0; if (geometry instanceof LineString) { const { coordinates } = geometry; return coordinates.slice(1).reduce((total, coord, i) => { const prev = coordinates[i]; return total + Point.distance(new Point(...prev), new Point(...coord)); }, 0); } if (geometry instanceof Polygon) { return geometry.rings.reduce( (perimeter, ring) => perimeter + Geometry.length(new LineString(ring)), 0, ); } throw new Error(`Unsupported geometry type: ${geometry.type}`); } /** * Find shortest connecting line between geometries * @param {Geometry} a - First geometry * @param {Geometry} b - Second geometry * @returns {LineString} Shortest connecting line */ static shortestLine(a, b) { const pointsA = Geometry.#geometryToPoints(a); const pointsB = Geometry.#geometryToPoints(b); let minDistance = Infinity; /** @type {[Point|null, Point|null]} */ let closestPair = [null, null]; for (const pointA of pointsA) { for (const pointB of pointsB) { const dist = Point.distance(pointA, pointB); if (dist < minDistance) { minDistance = dist; closestPair = [pointA, pointB]; } } } const [pointA, pointB] = closestPair; if (!pointA || !pointB) { throw new Error("Could not find shortest line between geometries"); } return new LineString([ [pointA.lng, pointA.lat], [pointB.lng, pointB.lat], ]); } /** * Check if geometries intersect * @param {Geometry} a - First geometry * @param {Geometry} b - Second geometry * @returns {boolean} */ static intersects(a, b) { if (!Geometry.#boundsIntersect(a, b)) return false; if (a instanceof Point) return Geometry.#pointIntersects(a, b); if (b instanceof Point) return Geometry.#pointIntersects(b, a); if (a instanceof LineString && b instanceof LineString) { return LineString.intersectsLine(a, b); } if (a instanceof Polygon) return Geometry.#polygonIntersects(a, b); if (b instanceof Polygon) return Geometry.#polygonIntersects(b, a); return Geometry.distance(a, b) < 0.1; } /** * Check if geometry A is within geometry B * @param {Geometry} a - Inner geometry * @param {Geometry} b - Outer geometry * @returns {boolean} */ static within(a, b) { const bA = Geometry.bounds(a); const bB = Geometry.bounds(b); // Quick bounds check if (bA.minX < bB.minX || bA.maxX > bB.maxX || bA.minY < bB.minY || bA.maxY > bB.maxY) { return false; } if (a instanceof Point) return Geometry.#pointWithin(a, b); if (a instanceof LineString && b instanceof Polygon) { return Geometry.#lineWithinPolygon(a, b); } if (a instanceof Polygon && b instanceof Polygon) { return Geometry.#polygonWithinPolygon(a, b); } return Geometry.#geometryToPoints(a).every((point) => Geometry.#pointWithin(point, b)); } /** * Check if geometry bounds intersect * @param {Geometry} a - First geometry * @param {Geometry} b - Second geometry * @returns {boolean} */ static #boundsIntersect(a, b) { const ba = Geometry.bounds(a); const bb = Geometry.bounds(b); return !(ba.maxX < bb.minX || ba.minX > bb.maxX || ba.maxY < bb.minY || ba.minY > bb.maxY); } /** * Check if point intersects geometry * @param {Point} point - Test point * @param {Geometry} geometry - Target geometry * @returns {boolean} */ static #pointIntersects(point, geometry) { if (geometry instanceof Point) return point.equals(geometry); if (geometry instanceof LineString) return Geometry.#pointOnLine(point, geometry); if (geometry instanceof Polygon) return Geometry.#pointInPolygon(point, geometry); return false; } /** * Check if point is within geometry * @param {Point} point - Test point * @param {Geometry} geometry - Container geometry * @returns {boolean} */ static #pointWithin(point, geometry) { if (geometry instanceof Point) return point.equals(geometry); if (geometry instanceof LineString) return Geometry.#pointOnLine(point, geometry); if (geometry instanceof Polygon) return Geometry.#pointInPolygon(point, geometry); return false; } /** * Check if point lies on line * @param {Point} point - Test point * @param {LineString} line - Target line * @param {number} [tolerance=0.1] - Tolerance in meters * @returns {boolean} */ static #pointOnLine(point, line, tolerance = 0.1) { const points = line.getPoints(); // Check vertices if (points.some((vertex) => point.equals(vertex, tolerance / 111320))) { return true; } // Check segments for (let i = 1; i < points.length; i++) { const p1 = points[i - 1]; const p2 = points[i]; const segmentLength = Point.distance(p1, p2); if (segmentLength === 0) continue; const d1 = Point.distance(point, p1); const d2 = Point.distance(point, p2); if (Math.abs(d1 + d2 - segmentLength) < tolerance) { return true; } } return false; } /** * Check if point is inside polygon using ray casting * @param {Point} point - Test point * @param {Polygon} polygon - Target polygon * @returns {boolean} */ static #pointInPolygon(point, polygon) { const [x, y] = [point.lng, point.lat]; const exterior = polygon.getExterior(); if (exterior.length < 3) return false; let inside = false; const n = exterior.length; for (let i = 0, j = n - 1; i < n; j = i++) { const [xi, yi] = exterior[i]; const [xj, yj] = exterior[j]; const intersect = yi > y !== yj > y && x < ((xj - xi) * (y - yi)) / (yj - yi) + xi; if (intersect) inside = !inside; } // Check holes if (inside) { for (const hole of polygon.getHoles()) { if (Geometry.#pointInRing(point, hole)) { return false; } } } return inside; } /** * Check if point is inside ring * @param {Point} point - Test point * @param {Coordinate[]} ring - Coordinate ring * @returns {boolean} */ static #pointInRing(point, ring) { const [x, y] = [point.lng, point.lat]; let inside = false; const n = ring.length; for (let i = 0, j = n - 1; i < n; j = i++) { const [xi, yi] = ring[i]; const [xj, yj] = ring[j]; const intersect = yi > y !== yj > y && x < ((xj - xi) * (y - yi)) / (yj - yi) + xi; if (intersect) inside = !inside; } return inside; } /** * Check if polygon intersects geometry * @param {Polygon} polygon - Test polygon * @param {Geometry} geometry - Target geometry * @returns {boolean} */ static #polygonIntersects(polygon, geometry) { // Check if any geometry points are inside polygon if ( Geometry.#geometryToPoints(geometry).some((point) => Geometry.#pointInPolygon(point, polygon)) ) { return true; } // Check intersections with exterior and holes const boundaries = [ new LineString(polygon.getExterior()), ...polygon.getHoles().map((hole) => new LineString(hole)), ]; return boundaries.some((boundary) => Geometry.intersects(boundary, geometry)); } /** * Check if line is within polygon * @param {LineString} line - Test line * @param {Polygon} polygon - Container polygon * @returns {boolean} */ static #lineWithinPolygon(line, polygon) { return ( line.getPoints().every((point) => Geometry.#pointInPolygon(point, polygon)) && !polygon.getHoles().some((hole) => Geometry.intersects(line, new LineString(hole))) ); } /** * Check if polygon A is within polygon B * @param {Polygon} polyA - Inner polygon * @param {Polygon} polyB - Outer polygon * @returns {boolean} */ static #polygonWithinPolygon(polyA, polyB) { return ( Geometry.#geometryToPoints(polyA).every((point) => Geometry.#pointInPolygon(point, polyB)) && !polyB.getHoles().some((hole) => Geometry.intersects(polyA, new LineString(hole))) ); } } /** Point geometry class */ export class Point extends Geometry { /** * @param {number} lng - Longitude * @param {number} lat - Latitude */ constructor(lng, lat) { super(); this.lng = lng; this.lat = lat; } /** @returns {Object} GeoJSON point */ toGeoJSON() { return { coordinates: [this.lng, this.lat], type: "Point", }; } /** @returns {Bounds} Point bounds */ bounds() { return new Bounds(this.lng, this.lat, this.lng, this.lat); } /** * Create Point from GeoJSON * @param {Object} geojson - GeoJSON point * @returns {Point} */ static fromGeoJSON(geojson) { const [lng, lat] = geojson.coordinates; return new Point(lng, lat); } /** * Calculate distance between points using Haversine formula * @param {Point} a - First point * @param {Point} b - Second point * @returns {number} Distance in meters */ static distance(a, b) { const φ1 = (a.lat * Math.PI) / 180; const φ2 = (b.lat * Math.PI) / 180; const Δφ = ((b.lat - a.lat) * Math.PI) / 180; const Δλ = ((b.lng - a.lng) * Math.PI) / 180; const aHav = Math.sin(Δφ / 2) ** 2 + Math.cos(φ1) * Math.cos(φ2) * Math.sin(Δλ / 2) ** 2; const c = 2 * Math.atan2(Math.sqrt(aHav), Math.sqrt(1 - aHav)); return R * c; } /** * Check if point equals another point * @param {Point} other - Comparison point * @param {number} [tolerance=TOLERANCE] - Tolerance in degrees * @returns {boolean} */ equals(other, tolerance = TOLERANCE) { return Math.abs(this.lng - other.lng) < tolerance && Math.abs(this.lat - other.lat) < tolerance; } } /** MultiLineString geometry class */ export class MultiLineString extends Geometry { /** * @param {Coordinate[][]} coordinates - Line coordinates */ constructor(coordinates) { super(); this.coordinates = coordinates; } /** * Create from GeoJSON * @param {Object} geojson - GeoJSON MultiLineString * @returns {MultiLineString} */ static fromGeoJSON(geojson) { return new MultiLineString(geojson.coordinates); } /** @returns {Object} GeoJSON representation */ toGeoJSON() { return { coordinates: this.coordinates, type: "MultiLineString", }; } /** @returns {Bounds} Geometry bounds */ bounds() { const allCoords = this.coordinates.flat(); return Geometry.calculateBoundsFromCoords(allCoords); } } /** LineString geometry class */ export class LineString extends Geometry { /** * @param {Coordinate[]} coordinates - Line coordinates */ constructor(coordinates) { super(); this.coordinates = coordinates; } /** * Check if two lines intersect * @param {LineString} line1 - First line * @param {LineString} line2 - Second line * @returns {boolean} */ static intersectsLine(line1, line2) { const coords1 = line1.coordinates; const coords2 = line2.coordinates; for (let i = 1; i < coords1.length; i++) { const [a1x, a1y] = coords1[i - 1]; const [a2x, a2y] = coords1[i]; for (let j = 1; j < coords2.length; j++) { const [b1x, b1y] = coords2[j - 1]; const [b2x, b2y] = coords2[j]; if (LineString.#segmentsIntersect(a1x, a1y, a2x, a2y, b1x, b1y, b2x, b2y)) { return true; } } } return false; } /** /** * Check if two segments intersect * @param {number} p0x * @param {number} p0y * @param {number} p1x * @param {number} p1y * @param {number} p2x * @param {number} p2y * @param {number} p3x * @param {number} p3y * @returns {boolean} */ static #segmentsIntersect(p0x, p0y, p1x, p1y, p2x, p2y, p3x, p3y) { const s1x = p1x - p0x; const s1y = p1y - p0y; const s2x = p3x - p2x; const s2y = p3y - p2y; const s = (-s1y * (p0x - p2x) + s1x * (p0y - p2y)) / (-s2x * s1y + s1x * s2y); const t = (s2x * (p0y - p2y) - s2y * (p0x - p2x)) / (-s2x * s1y + s1x * s2y); return s >= 0 && s <= 1 && t >= 0 && t <= 1; } /** * Interpolate point at distance along line * @param {LineString} line - Target line * @param {number} distance - Distance in meters * @param {Object} [options] - Interpolation options * @param {boolean} [options.normalized=false] - Treat distance as fraction * @returns {Point} Interpolated point */ static interpolatePoint(line, distance, options = {}) { const { normalized = false } = options; const coords = line.coordinates; if (coords.length < 2) { throw new Error("LineString must have at least 2 points"); } const targetDistance = normalized ? distance * Geometry.length(line) : distance; if (targetDistance < 0) return new Point(...coords[0]); let accumulatedDistance = 0; for (let i = 1; i < coords.length; i++) { const [lng1, lat1] = coords[i - 1]; const [lng2, lat2] = coords[i]; const segmentLength = Point.distance(new Point(lng1, lat1), new Point(lng2, lat2)); if (accumulatedDistance + segmentLength >= targetDistance) { const ratio = (targetDistance - accumulatedDistance) / segmentLength; const lng = lng1 + ratio * (lng2 - lng1); const lat = lat1 + ratio * (lat2 - lat1); return new Point(lng, lat); } accumulatedDistance += segmentLength; } return new Point(...coords[coords.length - 1]); } /** * Locate point on line and return distance * @param {LineString} line - Target line * @param {Point} point - Test point * @param {Object} [options] - Location options * @param {boolean} [options.normalized=false] - Return normalized distance * @returns {number} Distance in meters (or normalized) */ static locatePoint(line, point, options = {}) { const { normalized = false } = options; const coords = line.coordinates; if (coords.length < 2) { throw new Error("LineString must have at least 2 points"); } let minDistance = Infinity; let closestDistance = 0; let accumulatedDistance = 0; for (let i = 1; i < coords.length; i++) { const [lng1, lat1] = coords[i - 1]; const [lng2, lat2] = coords[i]; const segmentStart = new Point(lng1, lat1); const segmentEnd = new Point(lng2, lat2); const segmentLength = Point.distance(segmentStart, segmentEnd); if (segmentLength === 0) { accumulatedDistance += Point.distance(segmentStart, point); continue; } const u = ((point.lng - lng1) * (lng2 - lng1) + (point.lat - lat1) * (lat2 - lat1)) / ((lng2 - lng1) ** 2 + (lat2 - lat1) ** 2); const ratio = Math.max(0, Math.min(1, u)); const projectedLng = lng1 + ratio * (lng2 - lng1); const projectedLat = lat1 + ratio * (lat2 - lat1); const projectedPoint = new Point(projectedLng, projectedLat); const distToSegment = Point.distance(point, projectedPoint); if (distToSegment < minDistance) { minDistance = distToSegment; closestDistance = accumulatedDistance + ratio * segmentLength; } accumulatedDistance += segmentLength; } return normalized ? closestDistance / Geometry.length(line) : closestDistance; } /** @returns {Object} GeoJSON representation */ toGeoJSON() { return { coordinates: this.coordinates, type: "LineString", }; } /** @returns {Bounds} Line bounds */ bounds() { return Geometry.bounds(this); } /** * Create from GeoJSON * @param {Object} geojson - GeoJSON LineString * @returns {LineString} */ static fromGeoJSON(geojson) { return new LineString(geojson.coordinates); } /** * Get points as Point objects * @returns {Point[]} */ getPoints() { return this.coordinates.map(([lng, lat]) => new Point(lng, lat)); } } /** Polygon geometry class */ export class Polygon extends Geometry { /** * @param {Coordinate[][]} rings - Polygon rings (first is exterior, rest are holes) */ constructor(rings) { super(); this.rings = rings; } /** @returns {Object} GeoJSON representation */ toGeoJSON() { return { coordinates: this.rings, type: "Polygon", }; } /** @returns {Bounds} Polygon bounds */ bounds() { return Geometry.bounds(this); } /** * Calculate polygon centroid * @returns {Point} Centroid point */ centroid() { const exterior = this.getExterior(); let twiceArea = 0; let centroidLng = 0; let centroidLat = 0; for (let i = 0, j = exterior.length - 1; i < exterior.length; j = i++) { const [lng1, lat1] = exterior[j]; const [lng2, lat2] = exterior[i]; const cross = lng1 * lat2 - lng2 * lat1; twiceArea += cross; centroidLng += (lng1 + lng2) * cross; centroidLat += (lat1 + lat2) * cross; } const area = twiceArea * 0.5; if (Math.abs(area) < TOLERANCE) { // Degenerate polygon, return first point return new Point(...exterior[0]); } const factor = 1 / (6 * area); centroidLng *= factor; centroidLat *= factor; return new Point(centroidLng, centroidLat); } /** * Create from GeoJSON * @param {Object} geojson - GeoJSON Polygon * @returns {Polygon} */ static fromGeoJSON(geojson) { const rings = geojson.coordinates.map((ring /*: Coordinate[] */) => { if (ring.length < 1) return ring; const first = ring[0]; const last = ring[ring.length - 1]; if (first[0] !== last[0] || first[1] !== last[1]) { // Append first coordinate to the end if not closed return [...ring, first]; } return ring; }); return new Polygon(rings); } /** * Get exterior ring * @returns {Coordinate[]} */ getExterior() { return this.rings[0] ?? []; } /** * Get holes * @returns {Coordinate[][]} */ getHoles() { return this.rings.slice(1); } } /** GeoJSON Feature class */ export class Feature { /** * @param {Geometry} geometry - Feature geometry * @param {Object} [properties={}] - Feature properties * @param {string|number} [id] - Feature ID */ constructor(geometry, properties = {}, id = null) { this.type = "Feature"; this.geometry = geometry; this.properties = properties; this.id = id; } /** * Convert to GeoJSON object * @returns {Object} GeoJSON feature */ toGeoJSON() { return { geometry: this.geometry?.toGeoJSON() ?? null, properties: this.properties, type: "Feature", ...(this.id != null && { id: this.id }), }; } /** * Create from GeoJSON * @param {Object} geojson - GeoJSON feature * @returns {Feature|null} */ /** * @param {{geometry: object, properties?: object, id?: string|number}} geojson * @returns {Feature|null} */ static fromGeoJSON(geojson) { if (!geojson?.geometry) return null; const geometry = Geometry.fromGeoJSON(geojson.geometry); if (!geometry) return null; return new Feature(geometry, geojson.properties ?? {}, geojson.id ?? null); } } /** GeoJSON Feature Collection class */ export class Collection { /** * @param {Feature[]} [features=[]] - Feature array */ constructor(features = []) { this.type = "FeatureCollection"; this.features = features; } /** * Convert to GeoJSON object * @returns {Object} GeoJSON collection */ toGeoJSON() { return { features: this.features.map((feature) => feature.toGeoJSON()), type: "FeatureCollection", }; } /** * Create from GeoJSON * @param {Object} geojson - GeoJSON collection * @returns {Collection} */ /** * @param {{features?: any[]}} geojson * @returns {Collection} */ static fromGeoJSON(geojson) { const features = (geojson.features ?? []).map(Feature.fromGeoJSON).filter(Boolean); return new Collection(features); } /** * Get geometries by type * @param {string} type - Geometry type * @returns {Geometry[]} */ getGeometriesByType(type) { return this.features .filter((feature) => feature.geometry?.type === type) .map((feature) => feature.geometry); } } // Export constants export { R as EARTH_RADIUS, TOLERANCE }; // Convenience functions export const ser = (geometry) => geometry?.toGeoJSON() ?? null; export const de = (geojson) => Geometry.fromGeoJSON(geojson);