Add slim version of cheap-ruler impl

Author: zmiao

src/style-spec/expression/definitions/distance.js

@@ -15,9 +15,8 @@ import type {
     GeoJSONPolygon,
     GeoJSONMultiPolygon
 } from "@mapbox/geojson-types";
-import {classifyRings, updateBBox, boxWithinBox, pointWithinPolygon, segmentIntersectSegment} from '../../util/geometry_util.js';
-import CheapRuler from "cheap-ruler";
-import type {GeometryPoint} from "cheap-ruler";
+import {classifyRings, updateBBox, boxWithinBox, pointWithinPolygon, segmentIntersectSegment, CheapRuler} from '../../util/geometry_util.js';
+import type {GeometryPoint} from '../../util/geometry_util.js';
 import Point from "@mapbox/point-geometry";
 import Queue from "tinyqueue";
 

src/style-spec/util/geometry_util.js

@@ -3,7 +3,8 @@
 import quickselect from 'quickselect';
 
 import type Point from '@mapbox/point-geometry';
-import type {GeometryPoint} from "cheap-ruler";
+
+export type GeometryPoint = [number, number] | [number, number, number];
 /**
  * Returns the signed area for the polygon ring.  Postive areas are exterior rings and
  * have a clockwise winding.  Negative areas are interior rings and have a counter clockwise
@@ -145,3 +146,177 @@ export function segmentIntersectSegment(a: GeometryPoint, b: GeometryPoint, c: G
     if (twoSided(a, b, c, d) && twoSided(c, d, a, b)) return true;
     return false;
 }
+
+// normalize a degree value into [-180..180] range
+function wrap(deg) {
+    while (deg < -180) deg += 360;
+    while (deg > 180) deg -= 360;
+    return deg;
+}
+
+const factors = {
+    kilometers: 1,
+    miles: 1000 / 1609.344,
+    nauticalmiles: 1000 / 1852,
+    meters: 1000,
+    metres: 1000,
+    yards: 1000 / 0.9144,
+    feet: 1000 / 0.3048,
+    inches: 1000 / 0.0254
+};
+
+// Values that define WGS84 ellipsoid model of the Earth
+const RE = 6378.137; // equatorial radius
+const FE = 1 / 298.257223563; // flattening
+
+const E2 = FE * (2 - FE);
+const RAD = Math.PI / 180;
+
+/**
+ * A collection of very fast approximations to common geodesic measurements. Useful for performance-sensitive code that measures things on a city scale.
+ * This is the slim version of https://github.com/mapbox/cheap-ruler
+ * @param {number} lat latitude
+ * @param {string} [units='kilometers']
+ * @returns {CheapRuler}
+ * @example
+ * const ruler = cheapRuler(35.05, 'miles');
+ * //=ruler
+ */
+export class CheapRuler {
+    kx: number;
+    ky: number;
+    /**
+     * Creates a ruler instance for very fast approximations to common geodesic measurements around a certain latitude.
+     *
+     * @param {number} lat latitude
+     * @param {string} [units='kilometers']
+     * @returns {CheapRuler}
+     * @example
+     * const ruler = cheapRuler(35.05, 'miles');
+     * //=ruler
+     */
+    constructor(lat: number, units: string) {
+        if (lat === undefined) throw new Error('No latitude given.');
+        if (units && !factors[units]) throw new Error(`Unknown unit ${units}. Use one of: ${Object.keys(factors).join(', ')}`);
+
+        // Curvature formulas from https://en.wikipedia.org/wiki/Earth_radius#Meridional
+        const m = RAD * RE * (units ? factors[units] : 1);
+        const coslat = Math.cos(lat * RAD);
+        const w2 = 1 / (1 - E2 * (1 - coslat * coslat));
+        const w = Math.sqrt(w2);
+
+        // multipliers for converting longitude and latitude degrees into distance
+        this.kx = m * w * coslat;        // based on normal radius of curvature
+        this.ky = m * w * w2 * (1 - E2); // based on meridonal radius of curvature
+    }
+
+    /**
+     * Given two points of the form [longitude, latitude], returns the distance.
+     *
+     * @param {GeometryPoint} a point [longitude, latitude]
+     * @param {GeometryPoint} b point [longitude, latitude]
+     * @returns {number} distance
+     * @example
+     * const distance = ruler.distance([30.5, 50.5], [30.51, 50.49]);
+     * //=distance
+     */
+    distance(a: GeometryPoint, b: GeometryPoint): number {
+        const dx = wrap(a[0] - b[0]) * this.kx;
+        const dy = (a[1] - b[1]) * this.ky;
+        return Math.sqrt(dx * dx + dy * dy);
+    }
+
+    /**
+     * Returns the distance from a point `p` to a line segment `a` to `b`.
+     *
+     * @param {GeometryPoint} p point [longitude, latitude]
+     * @param {GeometryPoint} a segment point 1 [longitude, latitude]
+     * @param {GeometryPoint} b segment point 2 [longitude, latitude]
+     * @returns {number} distance
+     * @example
+     * const distance = ruler.pointToSegmentDistance([-67.04, 50.5], [-67.05, 50.57], [-67.03, 50.54]);
+     * //=distance
+     */
+    pointToSegmentDistance(p: GeometryPoint, a: GeometryPoint, b: GeometryPoint): number {
+        let [x, y] = a;
+        let dx = wrap(b[0] - x) * this.kx;
+        let dy = (b[1] - y) * this.ky;
+        let t = 0;
+
+        if (dx !== 0 || dy !== 0) {
+            t = (wrap(p[0] - x) * this.kx * dx + (p[1] - y) * this.ky * dy) / (dx * dx + dy * dy);
+
+            if (t > 1) {
+                x = b[0];
+                y = b[1];
+
+            } else if (t > 0) {
+                x += (dx / this.kx) * t;
+                y += (dy / this.ky) * t;
+            }
+        }
+
+        dx = wrap(p[0] - x) * this.kx;
+        dy = (p[1] - y) * this.ky;
+
+        return Math.sqrt(dx * dx + dy * dy);
+    }
+
+    /**
+     * Returns an object of the form {point, index, t}, where point is closest point on the line
+     * from the given point, index is the start index of the segment with the closest point,
+     * and t is a parameter from 0 to 1 that indicates where the closest point is on that segment.
+     *
+     * @param {Array<GeometryPoint>} line
+     * @param {GeometryPoint} p point [longitude, latitude]
+     * @returns {Object} {point, index, t}
+     * @example
+     * const point = ruler.pointOnLine(line, [-67.04, 50.5]).point;
+     * //=point
+     */
+    pointOnLine(line: Array<GeometryPoint>, p: GeometryPoint): Object {
+        let minDist = Infinity;
+        let minX = Infinity, minY = Infinity, minI = Infinity, minT  = Infinity;
+
+        for (let i = 0; i < line.length - 1; i++) {
+
+            let x = line[i][0];
+            let y = line[i][1];
+            let dx = wrap(line[i + 1][0] - x) * this.kx;
+            let dy = (line[i + 1][1] - y) * this.ky;
+            let t = 0;
+
+            if (dx !== 0 || dy !== 0) {
+                t = (wrap(p[0] - x) * this.kx * dx + (p[1] - y) * this.ky * dy) / (dx * dx + dy * dy);
+
+                if (t > 1) {
+                    x = line[i + 1][0];
+                    y = line[i + 1][1];
+
+                } else if (t > 0) {
+                    x += (dx / this.kx) * t;
+                    y += (dy / this.ky) * t;
+                }
+            }
+
+            dx = wrap(p[0] - x) * this.kx;
+            dy = (p[1] - y) * this.ky;
+
+            const sqDist = dx * dx + dy * dy;
+            if (sqDist < minDist) {
+                minDist = sqDist;
+                minX = x;
+                minY = y;
+                minI = i;
+                minT = t;
+            }
+        }
+
+        return {
+            point: [minX, minY],
+            index: minI,
+            t: Math.max(0, Math.min(1, minT))
+        };
+    }
+
+}