67 lines
2.5 KiB
TypeScript
67 lines
2.5 KiB
TypeScript
const R = 6371008.8;
|
||
const rad = (d: number) => (d * Math.PI) / 180;
|
||
|
||
export type LatLon = readonly [lat: number, lon: number];
|
||
|
||
export function haversine(a: LatLon, b: LatLon): number {
|
||
const dLat = rad(b[0] - a[0]);
|
||
const dLon = rad(b[1] - a[1]);
|
||
const s = Math.sin(dLat / 2) ** 2 + Math.cos(rad(a[0])) * Math.cos(rad(b[0])) * Math.sin(dLon / 2) ** 2;
|
||
return 2 * R * Math.asin(Math.min(1, Math.sqrt(s)));
|
||
}
|
||
|
||
/** Initial bearing a→b in degrees, 0 = north. */
|
||
export function bearing(a: LatLon, b: LatLon): number {
|
||
const y = Math.sin(rad(b[1] - a[1])) * Math.cos(rad(b[0]));
|
||
const x =
|
||
Math.cos(rad(a[0])) * Math.sin(rad(b[0])) -
|
||
Math.sin(rad(a[0])) * Math.cos(rad(b[0])) * Math.cos(rad(b[1] - a[1]));
|
||
return ((Math.atan2(y, x) * 180) / Math.PI + 360) % 360;
|
||
}
|
||
|
||
/** Distance in metres from point p to segment a–b (local equirectangular projection). */
|
||
export function distToSegment(p: LatLon, a: LatLon, b: LatLon): number {
|
||
const k = Math.cos(rad(p[0]));
|
||
const px = p[1] * k, py = p[0];
|
||
const ax = a[1] * k, ay = a[0];
|
||
const bx = b[1] * k, by = b[0];
|
||
const dx = bx - ax, dy = by - ay;
|
||
const len2 = dx * dx + dy * dy;
|
||
const t = len2 === 0 ? 0 : Math.max(0, Math.min(1, ((px - ax) * dx + (py - ay) * dy) / len2));
|
||
const cx = ax + t * dx, cy = ay + t * dy;
|
||
return Math.hypot(px - cx, py - cy) * (Math.PI / 180) * R;
|
||
}
|
||
|
||
/** Uniform grid for neighbourhood lookups; cell size in degrees of latitude. */
|
||
export class Grid<T> {
|
||
private cells = new Map<string, T[]>();
|
||
constructor(private cell = 0.0006) {}
|
||
private key(ix: number, iy: number) {
|
||
return `${ix}:${iy}`;
|
||
}
|
||
private ix(lon: number) {
|
||
return Math.floor(lon / (this.cell / Math.cos(rad(50))));
|
||
}
|
||
private iy(lat: number) {
|
||
return Math.floor(lat / this.cell);
|
||
}
|
||
add(lat: number, lon: number, item: T) {
|
||
const k = this.key(this.ix(lon), this.iy(lat));
|
||
const list = this.cells.get(k);
|
||
if (list) list.push(item);
|
||
else this.cells.set(k, [item]);
|
||
}
|
||
/** Items whose cell is within `radiusM` of the point (superset; caller refines). */
|
||
near(lat: number, lon: number, radiusM: number): T[] {
|
||
const dLat = radiusM / 111_320;
|
||
const x0 = this.ix(lon - dLat / Math.cos(rad(lat))), x1 = this.ix(lon + dLat / Math.cos(rad(lat)));
|
||
const y0 = this.iy(lat - dLat), y1 = this.iy(lat + dLat);
|
||
const out: T[] = [];
|
||
for (let x = x0; x <= x1; x++)
|
||
for (let y = y0; y <= y1; y++) {
|
||
const l = this.cells.get(this.key(x, y));
|
||
if (l) out.push(...l);
|
||
}
|
||
return out;
|
||
}
|
||
}
|