Simplify geometry

1 files changed, 143 insertions, 12 deletions

diff --git a/app/geometry.js b/app/geometry.js
index 62236e0..b9cac4f 100644
--- a/app/geometry.js
+++ b/app/geometry.js
@@ -9,7 +9,7 @@ const R = 6371000;
 const TOLERANCE = 1e-10;
 
 /**
- * Geographic bounds representation
+ * Geographic bounds representation for Bound caluclation typing
  */
 export class Bounds {
 	/**
@@ -186,6 +186,97 @@ export class Geometry {
 	}
 
 	/**
+	 * Calculate perpendicular distance from a point to a great-circle segment
+	 * @param {Point} point - Point to measure distance from
+	 * @param {Point} lineStart - Start point of the line segment
+	 * @param {Point} lineEnd - End point of the line segment
+	 * @returns {number} Perpendicular distance in meters
+	 */
+	static #perpendicularDistance(point, lineStart, lineEnd) {
+		// If start and end are the same, return distance to the point
+		if (lineStart.equals(lineEnd)) {
+			return Point.distance(point, lineStart);
+		}
+
+		// Convert coordinates to radians
+		const φ1 = (lineStart.lat * Math.PI) / 180;
+		const λ1 = (lineStart.lng * Math.PI) / 180;
+		const φ2 = (lineEnd.lat * Math.PI) / 180;
+		const λ2 = (lineEnd.lng * Math.PI) / 180;
+		const φ = (point.lat * Math.PI) / 180;
+		const λ = (point.lng * Math.PI) / 180;
+
+		// Calculate great-circle distance from start to end
+		const segmentLength = Point.distance(lineStart, lineEnd);
+		if (segmentLength === 0) {
+			return Point.distance(point, lineStart);
+		}
+
+		// Calculate bearing from lineStart to point and lineStart to lineEnd
+		const Δλ = λ2 - λ1;
+		const y = Math.sin(Δλ) * Math.cos(φ2);
+		const x = Math.cos(φ1) * Math.sin(φ2) - Math.sin(φ1) * Math.cos(φ2) * Math.cos(Δλ);
+		const bearing1 = Math.atan2(y, x);
+
+		const y2 = Math.sin(λ - λ1) * Math.cos(φ);
+		const x2 = Math.cos(φ1) * Math.sin(φ) - Math.sin(φ1) * Math.cos(φ) * Math.cos(λ - λ1);
+		const bearing2 = Math.atan2(y2, x2);
+
+		// Cross-track distance (approximation for small distances)
+		const crossTrack = Math.asin(
+			Math.sin(Point.distance(lineStart, point) / R) * Math.sin(bearing2 - bearing1),
+		);
+		const distance = Math.abs(R * crossTrack);
+
+		// Check if the closest point is within the segment
+		const alongTrack =
+			Math.acos(Math.cos(Point.distance(lineStart, point) / R) / Math.cos(distance / R)) * R;
+
+		if (alongTrack < 0 || alongTrack > segmentLength) {
+			// Return distance to closest endpoint
+			return Math.min(Point.distance(point, lineStart), Point.distance(point, lineEnd));
+		}
+
+		return distance;
+	}
+
+	/**
+	 * Simplify geometry using Ramer-Douglas-Peucker algorithm
+	 * @param {Coordinate[]} coords - Array of coordinates to simplify
+	 * @param {number} epsilon - Maximum distance tolerance in meters
+	 * @returns {Coordinate[]} Simplified coordinate array
+	 */
+	static simplifyCoords(coords, epsilon) {
+		if (coords.length <= 2) return coords;
+
+		let dmax = 0;
+		let index = 0;
+
+		// Find the point with the maximum perpendicular distance
+		const startPoint = new Point(coords[0][0], coords[0][1]);
+		const endPoint = new Point(coords[coords.length - 1][0], coords[coords.length - 1][1]);
+
+		for (let i = 1; i < coords.length - 1; i++) {
+			const point = new Point(coords[i][0], coords[i][1]);
+			const d = Geometry.#perpendicularDistance(point, startPoint, endPoint);
+			if (d > dmax) {
+				dmax = d;
+				index = i;
+			}
+		}
+
+		// If max distance is greater than epsilon, recursively simplify
+		if (dmax > epsilon) {
+			const recResults1 = Geometry.simplifyCoords(coords.slice(0, index + 1), epsilon);
+			const recResults2 = Geometry.simplifyCoords(coords.slice(index), epsilon);
+			// Combine results, excluding the duplicate point at index
+			return [...recResults1.slice(0, -1), ...recResults2];
+		} else {
+			return [coords[0], coords[coords.length - 1]];
+		}
+	}
+
+	/**
 	 * Calculate distance between geometries
 	 * @param {Geometry} a - First geometry
 	 * @param {Geometry} b - Second geometry
@@ -515,7 +606,6 @@ export class Geometry {
 	}
 }
 
-/** Point geometry class */
 export class Point extends Geometry {
 	/**
 	 * @param {number} lng - Longitude
@@ -579,7 +669,6 @@ export class Point extends Geometry {
 	}
 }
 
-/** MultiLineString geometry class */
 export class MultiLineString extends Geometry {
 	/**
 	 * @param {Coordinate[][]} coordinates - Line coordinates
@@ -606,6 +695,19 @@ export class MultiLineString extends Geometry {
 		};
 	}
 
+	/**
+	 * Simplify MultiLineString using Ramer-Douglas-Peucker algorithm
+	 * @param {number} epsilon - Maximum distance tolerance in meters
+	 * @returns {MultiLineString} Simplified MultiLineString
+	 */
+	simplify(epsilon) {
+		if (epsilon < 0) {
+			throw new Error("Epsilon must be non-negative");
+		}
+		const simplifiedCoords = this.coordinates.map((line) => Geometry.simplifyCoords(line, epsilon));
+		return new MultiLineString(simplifiedCoords);
+	}
+
 	/** @returns {Bounds} Geometry bounds */
 	bounds() {
 		const allCoords = this.coordinates.flat();
@@ -613,7 +715,6 @@ export class MultiLineString extends Geometry {
 	}
 }
 
-/** LineString geometry class */
 export class LineString extends Geometry {
 	/**
 	 * @param {Coordinate[]} coordinates - Line coordinates
@@ -624,6 +725,17 @@ export class LineString extends Geometry {
 	}
 
 	/**
+	 * @param {number} epsilon - Line coordinates
+	 */
+	simplify(epsilon) {
+		if (epsilon < 0) {
+			throw new Error("Epsilon must be non-negative");
+		}
+		const simplifiedCoords = Geometry.simplifyCoords(this.coordinates, epsilon);
+		return new LineString(simplifiedCoords);
+	}
+
+	/**
 	 * Check if two lines intersect
 	 * @param {LineString} line1 - First line
 	 * @param {LineString} line2 - Second line
@@ -798,7 +910,6 @@ export class LineString extends Geometry {
 	}
 }
 
-/** Polygon geometry class */
 export class Polygon extends Geometry {
 	/**
 	 * @param {Coordinate[][]} rings - Polygon rings (first is exterior, rest are holes)
@@ -857,6 +968,29 @@ export class Polygon extends Geometry {
 	}
 
 	/**
+	 * Simplify Polygon using Ramer-Douglas-Peucker algorithm
+	 * @param {number} epsilon - Maximum distance tolerance in meters
+	 * @returns {Polygon} Simplified Polygon
+	 */
+	simplify(epsilon) {
+		if (epsilon < 0) {
+			throw new Error("Epsilon must be non-negative");
+		}
+		const simplifiedRings = this.rings.map((ring) => {
+			const simplified = Geometry.simplifyCoords(ring, epsilon);
+			if (
+				simplified.length > 2 &&
+				(simplified[0][0] !== simplified[simplified.length - 1][0] ||
+					simplified[0][1] !== simplified[simplified.length - 1][1])
+			) {
+				simplified.push(simplified[0]);
+			}
+			return simplified;
+		});
+		return new Polygon(simplifiedRings);
+	}
+
+	/**
 	 * Create from GeoJSON
 	 * @param {Object} geojson - GeoJSON Polygon
 	 * @returns {Polygon}
@@ -892,14 +1026,13 @@ export class Polygon extends Geometry {
 	}
 }
 
-/** 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) {
+	constructor(geometry, properties = {}, id = "") {
 		this.type = "Feature";
 		this.geometry = geometry;
 		this.properties = properties;
@@ -934,11 +1067,10 @@ export class Feature {
 		const geometry = Geometry.fromGeoJSON(geojson.geometry);
 		if (!geometry) return null;
 
-		return new Feature(geometry, geojson.properties ?? {}, geojson.id ?? null);
+		return new Feature(geometry, geojson.properties ?? {}, geojson.id ?? "");
 	}
 }
 
-/** GeoJSON Feature Collection class */
 export class Collection {
 	/**
 	 * @param {Feature[]} [features=[]] - Feature array
@@ -965,12 +1097,11 @@ export class Collection {
 	 * @returns {Collection}
 	 */
 	/**
-	 * @param {{features?: any[]}} geojson
+	 * @param {{features: any[]}} geojson
 	 * @returns {Collection}
 	 */
 	static fromGeoJSON(geojson) {
-		const features = (geojson.features ?? []).map(Feature.fromGeoJSON).filter(Boolean);
-
+		const features = geojson.features ?? [];
 		return new Collection(features);
 	}