Add src/client/types/geoms.ts

src/client/types/geoms.ts

@@ -0,0 +1,177 @@
+// import L from "leaflet"
+
+const ellipsoid = {
+    a: 6378137,
+    b: 6356752.3142,
+    f: 1 / 298.257223563
+};
+
+function mod(n: number, p: number): number {
+    const r = n % p;
+    return r < 0 ? r + p : r;
+}
+
+function wrap(degrees: number, max = 360) {
+    if (-max <= degrees && degrees <= max) {
+        return degrees;
+    } else {
+        return mod(degrees + max, 2 * max) - max;
+    }
+}
+
+function dist(src, dst, itr = 100, mit = true): number {
+    const p1 = src,
+        p2 = dst;
+    const φ1 = toRadians(p1.lat),
+        λ1 = toRadians(p1.lng);
+    const φ2 = toRadians(p2.lat),
+        λ2 = toRadians(p2.lng);
+    const π = Math.PI;
+    const ε = Number.EPSILON;
+
+    // allow alternative ellipsoid to be specified
+    const { a, b, f } = ellipsoid;
+
+    const dL = λ2 - λ1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanφ.
+    const tanU1 = (1 - f) * Math.tan(φ1),
+        cosU1 = 1 / Math.sqrt(1 + tanU1 * tanU1),
+        sinU1 = tanU1 * cosU1;
+    const tanU2 = (1 - f) * Math.tan(φ2),
+        cosU2 = 1 / Math.sqrt(1 + tanU2 * tanU2),
+        sinU2 = tanU2 * cosU2;
+
+    const antipodal = Math.abs(dL) > π / 2 || Math.abs(φ2 - φ1) > π / 2;
+
+    let λ = dL,
+        sinλ: number | null = null,
+        cosλ: number | null = null; // λ = difference in longitude on an auxiliary sphere
+    let σ = antipodal ? π : 0,
+        sinσ = 0,
+        cosσ = antipodal ? -1 : 1,
+        sinSqσ: number | null = null; // σ = angular distance P₁ P₂ on the sphere
+    let cos2σₘ = 1; // σₘ = angular distance on the sphere from the equator to the midpoint of the line
+    let sinα: number | null = null,
+        cosSqα = 1; // α = azimuth of the geodesic at the equator
+    let C: number | null = null;
+
+    let λʹ: number | null = null,
+        iterations = 0;
+    do {
+        sinλ = Math.sin(λ);
+        cosλ = Math.cos(λ);
+        sinSqσ =
+            cosU2 * sinλ * (cosU2 * sinλ) +
+            (cosU1 * sinU2 - sinU1 * cosU2 * cosλ) * (cosU1 * sinU2 - sinU1 * cosU2 * cosλ);
+        if (Math.abs(sinSqσ) < ε) {
+            break; // co-incident/antipodal points (falls back on λ/σ = L)
+        }
+        sinσ = Math.sqrt(sinSqσ);
+        cosσ = sinU1 * sinU2 + cosU1 * cosU2 * cosλ;
+        σ = Math.atan2(sinσ, cosσ);
+        sinα = (cosU1 * cosU2 * sinλ) / sinσ;
+        cosSqα = 1 - sinα * sinα;
+        cos2σₘ = cosSqα !== 0 ? cosσ - (2 * sinU1 * sinU2) / cosSqα : 0; // on equatorial line cos²α = 0 (§6)
+        C = (f / 16) * cosSqα * (4 + f * (4 - 3 * cosSqα));
+        λʹ = λ;
+        λ = dL + (1 - C) * f * sinα * (σ + C * sinσ * (cos2σₘ + C * cosσ * (-1 + 2 * cos2σₘ * cos2σₘ)));
+        const iterationCheck = antipodal ? Math.abs(λ) - π : Math.abs(λ);
+        if (iterationCheck > π) {
+            throw new EvalError("λ > π");
+        }
+    } while (Math.abs(λ - λʹ) > 1e-12 && ++iterations < itr);
+
+    if (iterations >= itr) {
+        if (mit) {
+            return dist(
+                src,
+                new L.LatLng(dst.lat, dst.lng - 0.01),
+                itr,
+                mit
+            );
+        } else {
+            throw new EvalError(`Inverse vincenty formula failed to converge after ${itr} iterations 
+                    (start=${src.lat}/${src.lng}; dest=${dst.lat}/${dst.lng})`);
+        }
+    }
+    const uSq = (cosSqα * (a * a - b * b)) / (b * b);
+    const A = 1 + (uSq / 16384) * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
+    const B = (uSq / 1024) * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
+    const Δσ =
+        B *
+        sinσ *
+        (cos2σₘ +
+            (B / 4) *
+            (cosσ * (-1 + 2 * cos2σₘ * cos2σₘ) -
+                (B / 6) * cos2σₘ * (-3 + 4 * sinσ * sinσ) * (-3 + 4 * cos2σₘ * cos2σₘ)));
+
+    const s = b * A * (σ - Δσ); // s = length of the geodesic
+    return s
+
+}
+
+function pointDistance(src: L.LatLng, dst: L.LatLng): number {
+    return dist(
+        new L.LatLng(src.lat, wrap(src.lng, 180)),
+        new L.LatLng(dst.lat, wrap(dst.lng, 180))
+    );
+}
+
+function toRadians(degree: number): number {
+    return (degree * Math.PI) / 180;
+}
+
+function toDegrees(radians: number): number {
+    return (radians * 180) / Math.PI;
+}
+
+function midpoint(src: L.LatLng, dst: L.LatLng): L.LatLng {
+    // φm = atan2( sinφ1 + sinφ2, √( (cosφ1 + cosφ2⋅cosΔλ)² + cos²φ2⋅sin²Δλ ) )
+    // λm = λ1 + atan2(cosφ2⋅sinΔλ, cosφ1 + cosφ2⋅cosΔλ)
+    // midpoint is sum of vectors to two points: mathforum.org/library/drmath/view/51822.html
+
+    const φ1 = toRadians(src.lat);
+    const λ1 = toRadians(src.lng);
+    const φ2 = toRadians(dst.lat);
+    const Δλ = toRadians(dst.lng - src.lng);
+
+    // get cartesian coordinates for the two points
+    const A = { x: Math.cos(φ1), y: 0, z: Math.sin(φ1) }; // place point A on prime meridian y=0
+    const B = { x: Math.cos(φ2) * Math.cos(Δλ), y: Math.cos(φ2) * Math.sin(Δλ), z: Math.sin(φ2) };
+
+    // vector to midpoint is sum of vectors to two points (no need to normalise)
+    const C = { x: A.x + B.x, y: A.y + B.y, z: A.z + B.z };
+
+    const φm = Math.atan2(C.z, Math.sqrt(C.x * C.x + C.y * C.y));
+    const λm = λ1 + Math.atan2(C.y, C.x);
+
+    return new L.LatLng(toDegrees(φm), toDegrees(λm));
+}
+
+function recursiveMidPoint(src: L.LatLng, dst: L.LatLng, opt: { step?: number, dist?: number }): L.LatLng[] {
+    const geom: L.LatLng[] = [src, dst];
+    const mp = midpoint(src, dst);
+
+    if (opt.step != undefined) {
+        if (opt.step > 0) {
+            geom.splice(0, 1, ...recursiveMidPoint(src, mp, { step: opt.step - 1 }));
+            geom.splice(geom.length - 2, 2, ...recursiveMidPoint(mp, dst, { step: opt.step - 1 }));
+        } else {
+            geom.splice(1, 0, mp);
+        }
+    } else if (opt.dist != undefined) {
+        // console.log(src, dst, pointDistance(src, dst))
+        if (pointDistance(src, dst) > opt.dist) {
+            geom.splice(0, 1, ...recursiveMidPoint(src, mp, { dist: opt.dist }));
+            geom.splice(geom.length - 2, 2, ...recursiveMidPoint(mp, dst, { dist: opt.dist }));
+        } else {
+            geom.splice(1, 0, mp);
+        }
+    }
+    console.log(geom)
+    return geom;
+}
+
+export function getGeoLine(src: L.LatLng, dst: L.LatLng, opt: { step: number, dist: number }) {
+
+    return recursiveMidPoint(src, dst, opt)
+}