diff --git a/changelog.d/287-shared-great-circle-distance.md b/changelog.d/287-shared-great-circle-distance.md new file mode 100644 index 000000000..c2db9d54d --- /dev/null +++ b/changelog.d/287-shared-great-circle-distance.md @@ -0,0 +1,5 @@ +- The client's great-circle distances now come from one shared + `src/geoDistance.js` (`greatCircleKm`/`greatCircleMeters`, `{lat, lon}` + points). The Earth radius is standardized on 6371.0088 km, so distances shift + by about 1.4 ppm and region areas by about 2.8 ppm, and missing or + out-of-range points return `Infinity` (onionviolet, #287; Rajojha-1, #749). diff --git a/docs/CURRENT-STATE.md b/docs/CURRENT-STATE.md index 4c8f14f22..458c44c6e 100644 --- a/docs/CURRENT-STATE.md +++ b/docs/CURRENT-STATE.md @@ -1,5 +1,20 @@ # God's Eye View Current State +Shared distance checks use `src/geoDistance.js` with `{ lat, lon }` points and +the IUGG mean spherical Earth radius, 6371.0088 km. `greatCircleKm` is +canonical; `greatCircleMeters` is its exact metre conversion. Missing, +non-finite, and out-of-range points (latitude outside ±90, longitude outside +±180) return `Infinity`, so threshold and nearest-item consumers exclude invalid +positions. Older exported `haversineKm` entry points remain as compatibility +wrappers only. `src/data/approxDistance.js`, the equirectangular +viewer-proximity approximation, is deliberately separate, as are intentional +local angular, projection, ellipsoidal prefilter, and small-viewport +approximations. Four great-circle copies remain for a follow-up: +`src/data/adminBoundaries.js` and `src/tools/area.js` pass unwrapped longitudes +across the antimeridian, so they cannot use the strict ±180 check as-is, and +`src/sources/adsbRecords.js` and `server/providers/common/geo.js` keep their own +implementations. + Vessel sources can report healthy empty coverage using a current, complete snapshot with zero raw rows, a healthy transport and a positive last-message time. This clears obsolete contacts and settles startup without manufacturing diff --git a/scripts/package-boundaries.json b/scripts/package-boundaries.json index 40090ba19..43bc3d4de 100644 --- a/scripts/package-boundaries.json +++ b/scripts/package-boundaries.json @@ -301,6 +301,7 @@ "server/providers/regional/weather-effects.js", "server/providers/regional/weather.js", "src/data/regionalModel.js", + "src/geoDistance.js", "src/nominatimGeocode.js", "src/sources/httpBody.js", "src/i18n/marker.js", @@ -377,6 +378,7 @@ "src/search/nominatim.js", "src/nominatimGeocode.js", "src/data/regionalModel.js", + "src/geoDistance.js", "src/data/naturalEarthRegions.js", "src/data/bundledJson.js", "src/data/retryableLoad.js", @@ -915,6 +917,7 @@ "src/data/motionModel.js", "src/data/renderAltitude.js", "src/data/routePlausible.js", + "src/geoDistance.js", "src/data/trackedModelRegime.js", "src/data/trackingClickGesture.js", "src/layers/flights/controller.js", @@ -1009,6 +1012,7 @@ "src/data/motionModel.js", "src/data/renderAltitude.js", "src/data/routePlausible.js", + "src/geoDistance.js", "src/data/trackedModelRegime.js", "src/data/trackingClickGesture.js", "src/layers/aircraft/classification.js", @@ -1249,6 +1253,7 @@ "src/data/flowMatch.js", "src/data/tomtomTiles.js", "src/data/trafficBounds.js", + "src/geoDistance.js", "src/data/trafficFlowStyle.js", "src/data/trafficPresetStyle.js", "src/data/trafficQueue.js", @@ -1297,6 +1302,7 @@ "src/layers/bikeshare/ingestion.js", "src/layers/bikeshare/lifecycle.js", "src/layers/bikeshare/model.js", + "src/geoDistance.js", "src/sources/gbfsStations.js", "src/layers/bikeshare/policy.js", "src/layers/bikeshare/queries.js", @@ -1340,6 +1346,7 @@ "src/layers/cctv/ingestion.js", "src/layers/cctv/lifecycle.js", "src/layers/cctv/model.js", + "src/geoDistance.js", "src/layers/cctv/modelSpecs.js", "src/layers/cctv/navigation.js", "src/layers/cctv/policy.js", @@ -1527,6 +1534,7 @@ "exports": ["./geography/natural-earth", "./geography/neighborhoods"], "modules": [ "src/data/naturalEarthRegions.js", + "src/geoDistance.js", "src/data/neighborhoodPolygons.js", "src/data/bundledJson.js", "src/data/retryableLoad.js", @@ -1661,6 +1669,7 @@ "src/data/bhoteKoshiFloodPath.js", "src/data/bhoteKoshiShotPaths.js", "src/annotations/annotationEngine.js", + "src/geoDistance.js", "src/annotations/annotationResolver.js", "src/annotations/hybridAnnotationRenderer.js", "src/annotations/index.js", @@ -2513,6 +2522,7 @@ "src/data/bhoteKoshiIncidentPlaces.js", "src/data/bhoteKoshiFloodPath.js", "src/data/bhoteKoshiShotPaths.js", + "src/geoDistance.js", "src/annotations/annotationResolver.js", "src/annotations/resolver.js", "src/app/catalog.js", @@ -3075,6 +3085,7 @@ "src/data/approxDistance.js", "src/data/dataCredits.js", "src/data/geoid.js", + "src/geoDistance.js", "src/data/local_data/natural_earth/marine.json", "src/data/local_data/natural_earth/regions.json", "src/data/local_data/neighborhoods/san-francisco.json", @@ -3328,7 +3339,7 @@ }, "regional-source": { "exports": ["./sources/regional"], - "modules": ["src/data/regionalModel.js"], + "modules": ["src/data/regionalModel.js", "src/geoDistance.js"], "external": [] }, "voice-action-schemas": { @@ -3491,6 +3502,7 @@ "src/search/google.js", "src/search/geospatial.js", "src/data/regionalModel.js", + "src/geoDistance.js", "src/sources/httpBody.js", "src/i18n/marker.js", "src/sources/nominatim.js" @@ -3514,6 +3526,7 @@ "src/nominatimGeocode.js", "src/search/geospatial.js", "src/data/regionalModel.js", + "src/geoDistance.js", "src/sources/httpBody.js", "src/i18n/marker.js" ], @@ -3529,6 +3542,7 @@ "server/providers/transitHistory.js", "src/data/gtfsRealtime.js", "src/data/transitFeeds.js", + "src/geoDistance.js", "src/data/transitProxy.js", "src/sources/httpBody.js", "src/i18n/marker.js", @@ -3544,6 +3558,7 @@ "modules": [ "src/data/gtfsRealtime.js", "src/data/transitFeeds.js", + "src/geoDistance.js", "src/data/transitProxy.js", "src/sources/transitHistory.js" ], @@ -3565,6 +3580,7 @@ "src/data/geoid.js", "src/data/groundFloor.js", "src/data/gtfsRealtime.js", + "src/geoDistance.js", "src/data/iconOrientation.js", "src/data/inputOwnership.js", "src/data/labelArbiter.js", @@ -3622,6 +3638,7 @@ "src/data/iconOrientation.js", "src/data/inputOwnership.js", "src/data/transitFeeds.js", + "src/geoDistance.js", "src/data/transitIcons.js", "src/data/transitPresetStyle.js", "src/layers/transit/controls.js", @@ -3669,6 +3686,7 @@ "exports": ["./sources/transit-service"], "modules": [ "src/data/transitFeeds.js", + "src/geoDistance.js", "src/data/transitProxy.js", "src/data/gtfsRealtime.js", "src/sources/httpBody.js", @@ -3793,7 +3811,8 @@ "src/view/index.js", "src/tools/views.js", "src/tools/globePanel.js", - "src/tools/queries/panelRequest.js" + "src/tools/queries/panelRequest.js", + "src/geoDistance.js" ], "external": ["satellite.js"] }, diff --git a/src/annotations/annotationEngine.js b/src/annotations/annotationEngine.js index 645aa3f3d..b0622579d 100644 --- a/src/annotations/annotationEngine.js +++ b/src/annotations/annotationEngine.js @@ -1,6 +1,7 @@ import { isUnavailableCapability } from '../sources/capability.js'; import { defaultGeospatial } from '../search/defaults.js'; import * as Cesium from 'cesium'; +import { greatCircleMeters } from '../geoDistance.js'; import { holdContinuousRender, releaseContinuousRender, @@ -1567,16 +1568,7 @@ async function fetchRoute(coordPairs, mode, signal, service) { } function greatCircleM(a, b) { - const R = 6371000; - const toRad = (deg) => (deg * Math.PI) / 180; - const dLat = toRad(b.lat - a.lat); - const dLon = toRad(b.lon - a.lon); - const la1 = toRad(a.lat); - const la2 = toRad(b.lat); - const h = - Math.sin(dLat / 2) ** 2 + - Math.cos(la1) * Math.cos(la2) * Math.sin(dLon / 2) ** 2; - return 2 * R * Math.asin(Math.min(1, Math.sqrt(h))); + return greatCircleMeters(a, b); } function formatDistance(m) { diff --git a/src/annotations/drawMode.js b/src/annotations/drawMode.js index 34f08e8f1..02a501010 100644 --- a/src/annotations/drawMode.js +++ b/src/annotations/drawMode.js @@ -20,6 +20,8 @@ import { t, tn } from '../i18n/index.js'; * `drawTool.js`. */ +import { greatCircleMeters } from '../geoDistance.js'; + export const DRAW_SHAPES = Object.freeze(['area', 'line', 'pin']); export const MIN_VERTICES = Object.freeze({ area: 3, line: 2, pin: 1 }); /** @@ -49,14 +51,7 @@ export function createDrawSession(shape = 'area') { /** Great-circle distance in metres between two {lon, lat} points. */ export function greatCircleM(a, b) { - const R = 6371000; - const toRad = (d) => (d * Math.PI) / 180; - const dLat = toRad(b.lat - a.lat); - const dLon = toRad(b.lon - a.lon); - const h = - Math.sin(dLat / 2) ** 2 + - Math.cos(toRad(a.lat)) * Math.cos(toRad(b.lat)) * Math.sin(dLon / 2) ** 2; - return 2 * R * Math.asin(Math.sqrt(Math.min(1, h))); + return greatCircleMeters(a, b); } /** diff --git a/src/cockpitCloudEffects.js b/src/cockpitCloudEffects.js index 053217810..ccbf31831 100644 --- a/src/cockpitCloudEffects.js +++ b/src/cockpitCloudEffects.js @@ -1,5 +1,6 @@ import { applicationServices } from './services/application.js'; import * as Cesium from 'cesium'; +import { greatCircleMeters } from './geoDistance.js'; import { deriveWeatherEffectProfile, weatherAltitudeFactors, @@ -127,23 +128,9 @@ function compileShader(gl, type, source) { } function greatCircleM(a, b) { - if ( - ![a?.latitude, a?.longitude, b?.latitude, b?.longitude].every( - Number.isFinite, - ) - ) - return Infinity; - const latitudeA = Cesium.Math.toRadians(a.latitude); - const latitudeB = Cesium.Math.toRadians(b.latitude); - const latitudeDelta = Cesium.Math.toRadians(b.latitude - a.latitude); - const longitudeDelta = Cesium.Math.toRadians(b.longitude - a.longitude); - const haversine = - Math.sin(latitudeDelta / 2) ** 2 + - Math.cos(latitudeA) * - Math.cos(latitudeB) * - Math.sin(longitudeDelta / 2) ** 2; - return ( - 6371000 * 2 * Math.atan2(Math.sqrt(haversine), Math.sqrt(1 - haversine)) + return greatCircleMeters( + { lat: a?.latitude, lon: a?.longitude }, + { lat: b?.latitude, lon: b?.longitude }, ); } diff --git a/src/data/analystEngine.js b/src/data/analystEngine.js index 697469dba..6c3361aff 100644 --- a/src/data/analystEngine.js +++ b/src/data/analystEngine.js @@ -27,6 +27,7 @@ import { feedProvenanceEnvelope } from './layerSnapshot.js'; import { pointInRing } from './naturalEarthRegions.js'; import { VOICE_LAYER_MANIFEST, voiceLayer } from '../voice/layerManifest.js'; +import { greatCircleKm } from '../geoDistance.js'; /** * Layers the engine understands, with the fields queries may reference, @@ -98,17 +99,9 @@ function toEpochMs(value) { return Date.parse(text); } -const EARTH_R_KM = 6371; - -/** Great-circle distance in km. */ +/** @deprecated Use geoDistance.greatCircleKm with `{ lat, lon }` points. */ export function haversineKm(lat1, lon1, lat2, lon2) { - const d2r = Math.PI / 180; - const dLat = (lat2 - lat1) * d2r; - const dLon = (lon2 - lon1) * d2r; - const a = - Math.sin(dLat / 2) ** 2 + - Math.cos(lat1 * d2r) * Math.cos(lat2 * d2r) * Math.sin(dLon / 2) ** 2; - return 2 * EARTH_R_KM * Math.asin(Math.min(1, Math.sqrt(a))); + return greatCircleKm({ lat: lat1, lon: lon1 }, { lat: lat2, lon: lon2 }); } /** diff --git a/src/data/naturalEarthRegions.js b/src/data/naturalEarthRegions.js index 0c7af42d0..c40e180ac 100644 --- a/src/data/naturalEarthRegions.js +++ b/src/data/naturalEarthRegions.js @@ -18,8 +18,8 @@ import { loadBundledJson } from './bundledJson.js'; import { createRetryableLoader } from './retryableLoad.js'; +import { EARTH_RADIUS_KM, greatCircleKm } from '../geoDistance.js'; -const EARTH_RADIUS_KM = 6371; const toRad = (d) => (d * Math.PI) / 180; /** Spherical-excess ring area (km²) — same family as turf/geojson-area. */ @@ -36,15 +36,6 @@ function ringAreaKm2(ring) { return Math.abs((sum * EARTH_RADIUS_KM * EARTH_RADIUS_KM) / 2); } -function haversineKm(lon1, lat1, lon2, lat2) { - const dLat = toRad(lat2 - lat1); - const dLon = toRad(lon2 - lon1); - const h = - Math.sin(dLat / 2) ** 2 + - Math.cos(toRad(lat1)) * Math.cos(toRad(lat2)) * Math.sin(dLon / 2) ** 2; - return 2 * EARTH_RADIUS_KM * Math.asin(Math.min(1, Math.sqrt(h))); -} - /** * Normalize a query/name for matching: lowercase, strip diacritics and * punctuation, collapse whitespace, strip a leading "the ". @@ -157,7 +148,10 @@ function buildEntries(pack, kind) { polygons, areaKm2, bbox: [minLon, minLat, maxLon, maxLat], - bboxDiagonalKm: haversineKm(minLon, minLat, maxLon, maxLat), + bboxDiagonalKm: greatCircleKm( + { lat: minLat, lon: minLon }, + { lat: maxLat, lon: maxLon }, + ), }); } return out; diff --git a/src/data/naturalEarthRegions.test.mjs b/src/data/naturalEarthRegions.test.mjs index 36c1946ad..5dd78adea 100644 --- a/src/data/naturalEarthRegions.test.mjs +++ b/src/data/naturalEarthRegions.test.mjs @@ -3,6 +3,7 @@ import { test } from 'node:test'; import assert from 'node:assert/strict'; import { statSync, readFileSync } from 'node:fs'; import { findNaturalRegion, listRegions, lookupNaturalRegionOutline, pointInRing } from './naturalEarthRegions.js'; +import { greatCircleKm } from '../geoDistance.js'; const PACK_DIR = new URL('./local_data/natural_earth/', import.meta.url); @@ -19,6 +20,16 @@ test('marquee ranges resolve with sane areas (owner acceptance: Alps + Rockies)' assert.ok(rockies.areaKm2 >= 700_000 && rockies.areaKm2 < 1_500_000, `Rockies area sane (got ${Math.round(rockies.areaKm2)} km²)`); assert.ok(rockies.bboxDiagonalKm > 2000, 'Rockies span thousands of km'); + const coordinates = rockies.polygons.flat(); + const minLon = Math.min(...coordinates.map(([lon]) => lon)); + const minLat = Math.min(...coordinates.map(([, lat]) => lat)); + const maxLon = Math.max(...coordinates.map(([lon]) => lon)); + const maxLat = Math.max(...coordinates.map(([, lat]) => lat)); + assert.equal( + rockies.bboxDiagonalKm, + greatCircleKm({ lat: minLat, lon: minLon }, { lat: maxLat, lon: maxLon }), + 'Natural Earth bbox conversion maps [lon, lat] data into { lat, lon } points', + ); }); test('aliases and articles: "the Alps", "Rockies", "Sahara Desert", "Himalaya"', async () => { diff --git a/src/data/regionalModel.js b/src/data/regionalModel.js index fca473994..46f162d6f 100644 --- a/src/data/regionalModel.js +++ b/src/data/regionalModel.js @@ -1,3 +1,5 @@ +import { greatCircleMeters } from '../geoDistance.js'; + const MAX_ARTICLES = 5; function cleanText(value, maxLength = 180) { @@ -153,19 +155,8 @@ export function weatherCodeLabel(code) { /** Great-circle distance used to avoid refetching a regional brief every animation frame. */ export function regionalDistanceM(from, to) { - if ( - ![from?.latitude, from?.longitude, to?.latitude, to?.longitude].every( - Number.isFinite, - ) - ) { - return Infinity; - } - const phi1 = (from.latitude * Math.PI) / 180; - const phi2 = (to.latitude * Math.PI) / 180; - const deltaPhi = ((to.latitude - from.latitude) * Math.PI) / 180; - const deltaLambda = ((to.longitude - from.longitude) * Math.PI) / 180; - const a = - Math.sin(deltaPhi / 2) ** 2 + - Math.cos(phi1) * Math.cos(phi2) * Math.sin(deltaLambda / 2) ** 2; - return 6371000 * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); + return greatCircleMeters( + { lat: from?.latitude, lon: from?.longitude }, + { lat: to?.latitude, lon: to?.longitude }, + ); } diff --git a/src/data/routePlausible.js b/src/data/routePlausible.js index a8c0e9259..238c2e6cb 100644 --- a/src/data/routePlausible.js +++ b/src/data/routePlausible.js @@ -10,19 +10,9 @@ * low → its destination should be nearby). */ -const D2R = Math.PI / 180; -const R_KM = 6371; +import { EARTH_RADIUS_KM, greatCircleKm, isGeoPoint } from '../geoDistance.js'; -/** Haversine great-circle distance in km. */ -export function greatCircleKm(lat1, lon1, lat2, lon2) { - const p1 = lat1 * D2R; - const p2 = lat2 * D2R; - const dp = (lat2 - lat1) * D2R; - const dl = (lon2 - lon1) * D2R; - const a = - Math.sin(dp / 2) ** 2 + Math.cos(p1) * Math.cos(p2) * Math.sin(dl / 2) ** 2; - return R_KM * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); -} +const D2R = Math.PI / 180; function bearingRad(lat1, lon1, lat2, lon2) { const p1 = lat1 * D2R; @@ -36,10 +26,11 @@ function bearingRad(lat1, lon1, lat2, lon2) { /** Signed cross-track distance (km) of point from the great circle p1→p2. */ export function crossTrackKm(lat, lon, lat1, lon1, lat2, lon2) { - const d13 = greatCircleKm(lat1, lon1, lat, lon) / R_KM; + const d13 = + greatCircleKm({ lat: lat1, lon: lon1 }, { lat, lon }) / EARTH_RADIUS_KM; const b13 = bearingRad(lat1, lon1, lat, lon); const b12 = bearingRad(lat1, lon1, lat2, lon2); - return Math.asin(Math.sin(d13) * Math.sin(b13 - b12)) * R_KM; + return Math.asin(Math.sin(d13) * Math.sin(b13 - b12)) * EARTH_RADIUS_KM; } const NEAR_ENDPOINT_KM = 130; @@ -65,14 +56,14 @@ export function routePlausible({ origin = null, destination = null, }) { - const haveO = Number.isFinite(origin?.lat) && Number.isFinite(origin?.lon); - const haveD = - Number.isFinite(destination?.lat) && Number.isFinite(destination?.lon); + const position = { lat: latDeg, lon: lonDeg }; + if (!isGeoPoint(position)) return true; + const haveO = isGeoPoint(origin); + const haveD = isGeoPoint(destination); if (!haveO && !haveD) return true; // no coordinates — cannot judge, do not hide // (a) Geographic consistency: near an endpoint, or roughly on the path. - const near = (pt) => - greatCircleKm(latDeg, lonDeg, pt.lat, pt.lon) < NEAR_ENDPOINT_KM; + const near = (point) => greatCircleKm(position, point) < NEAR_ENDPOINT_KM; let geomOk = (haveO && near(origin)) || (haveD && near(destination)); if (!geomOk && haveO && haveD) { geomOk = @@ -99,17 +90,10 @@ export function routePlausible({ Math.abs(verticalRateMps) > VERT_TREND_MPS ) { if (verticalRateMps > 0) { - if ( - haveO && - greatCircleKm(latDeg, lonDeg, origin.lat, origin.lon) > LOCAL_AIRPORT_KM - ) + if (haveO && greatCircleKm(position, origin) > LOCAL_AIRPORT_KM) return false; // departing — origin should be local } else { - if ( - haveD && - greatCircleKm(latDeg, lonDeg, destination.lat, destination.lon) > - LOCAL_AIRPORT_KM - ) + if (haveD && greatCircleKm(position, destination) > LOCAL_AIRPORT_KM) return false; // arriving — destination should be local } } diff --git a/src/data/routePlausible.test.mjs b/src/data/routePlausible.test.mjs index dab932c6e..aa59a0c18 100644 --- a/src/data/routePlausible.test.mjs +++ b/src/data/routePlausible.test.mjs @@ -1,14 +1,15 @@ // src/data/routePlausible.test.mjs import { test } from 'node:test'; import assert from 'node:assert/strict'; -import { routePlausible, greatCircleKm } from './routePlausible.js'; +import { greatCircleKm } from '../geoDistance.js'; +import { routePlausible } from './routePlausible.js'; const SFO = { lat: 37.6188, lon: -122.3754 }; const LAX = { lat: 33.9416, lon: -118.4085 }; const JFK = { lat: 40.6413, lon: -73.7781 }; test('greatCircleKm sanity: SFO→LAX ≈ 543 km', () => { - const d = greatCircleKm(SFO.lat, SFO.lon, LAX.lat, LAX.lon); + const d = greatCircleKm(SFO, LAX); assert.ok(Math.abs(d - 543) < 15, `got ${d}`); }); @@ -41,4 +42,8 @@ test('missing data never hides a route (cannot judge → allow)', () => { latDeg: 51.5, lonDeg: -0.12, origin: { lat: null, lon: null }, destination: null, }), true); + assert.equal(routePlausible({ + latDeg: 95, lonDeg: 0, + origin: SFO, destination: LAX, + }), true); }); diff --git a/src/data/trafficBounds.js b/src/data/trafficBounds.js index ff5142133..5919f26f1 100644 --- a/src/data/trafficBounds.js +++ b/src/data/trafficBounds.js @@ -13,31 +13,11 @@ * @module data/trafficBounds */ -/** @const {number} Mean Earth radius in km (spherical approximation). */ -const EARTH_RADIUS_KM = 6371; +import { EARTH_RADIUS_KM, greatCircleKm, isGeoPoint } from '../geoDistance.js'; const toRad = (deg) => (deg * Math.PI) / 180; const toDeg = (rad) => (rad * 180) / Math.PI; -/** - * Great-circle distance between two lat/lon points (haversine). - * - * @param {number} lat1 - First point latitude (degrees). - * @param {number} lon1 - First point longitude (degrees). - * @param {number} lat2 - Second point latitude (degrees). - * @param {number} lon2 - Second point longitude (degrees). - * @returns {number} Distance in kilometres. - */ -export function greatCircleKm(lat1, lon1, lat2, lon2) { - const p1 = toRad(lat1); - const p2 = toRad(lat2); - const dp = toRad(lat2 - lat1); - const dl = toRad(lon2 - lon1); - const a = - Math.sin(dp / 2) ** 2 + Math.cos(p1) * Math.cos(p2) * Math.sin(dl / 2) ** 2; - return 2 * EARTH_RADIUS_KM * Math.asin(Math.min(1, Math.sqrt(a))); -} - /** * Initial bearing (radians) from point 1 toward point 2 along the great circle. * @@ -111,10 +91,13 @@ export function deriveFetchCenter({ hitLon, maxPullKm = 12, }) { - if (!Number.isFinite(hitLat) || !Number.isFinite(hitLon)) { + if (!isGeoPoint({ lat: hitLat, lon: hitLon })) { return { lat: nadirLat, lon: nadirLon, source: 'nadir' }; } - const distKm = greatCircleKm(nadirLat, nadirLon, hitLat, hitLon); + const distKm = greatCircleKm( + { lat: nadirLat, lon: nadirLon }, + { lat: hitLat, lon: hitLon }, + ); if (distKm <= maxPullKm) { return { lat: hitLat, lon: hitLon, source: 'hit' }; } diff --git a/src/data/trafficBounds.test.mjs b/src/data/trafficBounds.test.mjs index 0ff38628d..bbc60cbfe 100644 --- a/src/data/trafficBounds.test.mjs +++ b/src/data/trafficBounds.test.mjs @@ -1,8 +1,8 @@ // src/data/trafficBounds.test.mjs import { test } from 'node:test'; import assert from 'node:assert/strict'; +import { greatCircleKm } from '../geoDistance.js'; import { - greatCircleKm, deriveFetchCenter, clampBoundsAroundCenter, } from './trafficBounds.js'; @@ -11,7 +11,7 @@ import { const NADIR = { lat: 30.2672, lon: -97.7431 }; test('greatCircleKm sanity: Austin -> ~5 km north', () => { - const d = greatCircleKm(NADIR.lat, NADIR.lon, NADIR.lat + 0.045, NADIR.lon); + const d = greatCircleKm(NADIR, { lat: NADIR.lat + 0.045, lon: NADIR.lon }); assert.ok(Math.abs(d - 5.0) < 0.1, `got ${d}`); }); @@ -41,14 +41,14 @@ test('oblique look within 12 km: uses the hit point verbatim', () => { test('horizon gaze: far hit is pulled back to 12 km along the bearing', () => { // Hit ~100 km due east of nadir. const hit = { lat: NADIR.lat, lon: NADIR.lon + 1.041 }; - assert.ok(greatCircleKm(NADIR.lat, NADIR.lon, hit.lat, hit.lon) > 90, 'precondition: far hit'); + assert.ok(greatCircleKm(NADIR, hit) > 90, 'precondition: far hit'); const c = deriveFetchCenter({ nadirLat: NADIR.lat, nadirLon: NADIR.lon, hitLat: hit.lat, hitLon: hit.lon, maxPullKm: 12, }); assert.equal(c.source, 'pulled'); - const d = greatCircleKm(NADIR.lat, NADIR.lon, c.lat, c.lon); + const d = greatCircleKm(NADIR, c); assert.ok(Math.abs(d - 12) < 0.05, `pulled distance ${d} km, expected ~12`); // Due-east bearing: latitude stays ~constant, longitude moves east but well // short of the hit. @@ -64,12 +64,12 @@ test('pull cap is honored for other maxPullKm values', () => { maxPullKm: 5, }); assert.equal(c.source, 'pulled'); - const d = greatCircleKm(NADIR.lat, NADIR.lon, c.lat, c.lon); + const d = greatCircleKm(NADIR, c); assert.ok(Math.abs(d - 5) < 0.05, `pulled distance ${d} km, expected ~5`); }); -test('pickEllipsoid failure (non-finite hit): falls back to the camera nadir', () => { - for (const [hitLat, hitLon] of [[NaN, NaN], [undefined, undefined], [30.3, undefined]]) { +test('invalid ground hit falls back to the camera nadir', () => { + for (const [hitLat, hitLon] of [[NaN, NaN], [undefined, undefined], [30.3, undefined], [91, -97], [30, 181]]) { const c = deriveFetchCenter({ nadirLat: NADIR.lat, nadirLon: NADIR.lon, hitLat, hitLon, diff --git a/src/data/transitFeeds.js b/src/data/transitFeeds.js index dfd3d6b18..b8476a6ee 100644 --- a/src/data/transitFeeds.js +++ b/src/data/transitFeeds.js @@ -1,3 +1,5 @@ +import { greatCircleKm } from '../geoDistance.js'; + /** * @module transitFeeds * @description Registry of keyless, openly licensed GTFS-Realtime @@ -773,7 +775,8 @@ export function getRegisteredTransitFeed(id) { } /** - * Great-circle distance in km. + * Great-circle distance in km. Retained for compatibility with existing + * transit consumers; new shared spatial code should use `greatCircleKm`. * @param {number} aLat * @param {number} aLon * @param {number} bLat @@ -781,13 +784,7 @@ export function getRegisteredTransitFeed(id) { * @returns {number} */ export function haversineKm(aLat, aLon, bLat, bLon) { - const toRad = (deg) => (deg * Math.PI) / 180; - const dLat = toRad(bLat - aLat); - const dLon = toRad(bLon - aLon); - const h = - Math.sin(dLat / 2) ** 2 + - Math.cos(toRad(aLat)) * Math.cos(toRad(bLat)) * Math.sin(dLon / 2) ** 2; - return 2 * 6371 * Math.asin(Math.min(1, Math.sqrt(h))); + return greatCircleKm({ lat: aLat, lon: aLon }, { lat: bLat, lon: bLon }); } /** diff --git a/src/geoDistance.js b/src/geoDistance.js new file mode 100644 index 000000000..e71b72bb3 --- /dev/null +++ b/src/geoDistance.js @@ -0,0 +1,70 @@ +/** + * Canonical geographic point used by shared spatial utilities. + * + * @typedef {{lat: number, lon: number}} GeoPoint + */ + +/** IUGG mean spherical Earth radius, in kilometres. */ +export const EARTH_RADIUS_KM = 6371.0088; + +const DEGREES_TO_RADIANS = Math.PI / 180; + +/** + * Whether a value is a canonical, finite WGS84-style point in degrees. + * Longitudes are not implicitly wrapped because an out-of-range value usually + * indicates a broken coordinate contract at the caller. + * + * @param {*} point + * @returns {point is GeoPoint} + */ +export function isGeoPoint(point) { + return ( + Number.isFinite(point?.lat) && + point.lat >= -90 && + point.lat <= 90 && + Number.isFinite(point?.lon) && + point.lon >= -180 && + point.lon <= 180 + ); +} + +/** + * Great-circle distance in kilometres between two `{ lat, lon }` points. + * + * `{ lat, lon }` is the canonical shared coordinate representation. Invalid, + * missing, non-finite, or out-of-range coordinates return `Infinity`, which + * keeps invalid points out of threshold and nearest-neighbour selections. + * + * @param {GeoPoint} from + * @param {GeoPoint} to + * @returns {number} + */ +export function greatCircleKm(from, to) { + if (!isGeoPoint(from) || !isGeoPoint(to)) return Number.POSITIVE_INFINITY; + + const latitude1 = from.lat * DEGREES_TO_RADIANS; + const latitude2 = to.lat * DEGREES_TO_RADIANS; + const latitudeDelta = (to.lat - from.lat) * DEGREES_TO_RADIANS; + const longitudeDelta = (to.lon - from.lon) * DEGREES_TO_RADIANS; + const haversine = + Math.sin(latitudeDelta / 2) ** 2 + + Math.cos(latitude1) * + Math.cos(latitude2) * + Math.sin(longitudeDelta / 2) ** 2; + const bounded = Math.min(1, Math.max(0, haversine)); + return ( + 2 * EARTH_RADIUS_KM * Math.atan2(Math.sqrt(bounded), Math.sqrt(1 - bounded)) + ); +} + +/** + * Great-circle distance in metres, with the same point and invalid-input + * contract as {@link greatCircleKm}. + * + * @param {GeoPoint} from + * @param {GeoPoint} to + * @returns {number} + */ +export function greatCircleMeters(from, to) { + return greatCircleKm(from, to) * 1000; +} diff --git a/src/geoDistance.test.mjs b/src/geoDistance.test.mjs new file mode 100644 index 000000000..67d6c32fc --- /dev/null +++ b/src/geoDistance.test.mjs @@ -0,0 +1,78 @@ +import test from 'node:test'; +import assert from 'node:assert/strict'; +import { + EARTH_RADIUS_KM, + greatCircleKm, + greatCircleMeters, + isGeoPoint, +} from './geoDistance.js'; + +const AUSTIN = Object.freeze({ lat: 30.2672, lon: -97.7431 }); +const HOUSTON = Object.freeze({ lat: 29.7604, lon: -95.3698 }); + +test('great-circle identity and symmetry', () => { + assert.equal(greatCircleKm(AUSTIN, AUSTIN), 0); + assert.equal(greatCircleKm(AUSTIN, HOUSTON), greatCircleKm(HOUSTON, AUSTIN)); +}); + +test('great-circle one-degree calibration and canonical units', () => { + const kilometres = greatCircleKm({ lat: 0, lon: 0 }, { lat: 1, lon: 0 }); + assert.ok(Math.abs(kilometres - 111.195) < 0.001); + assert.equal( + greatCircleMeters({ lat: 0, lon: 0 }, { lat: 1, lon: 0 }), + kilometres * 1000, + ); +}); + +test('great-circle handles dateline crossing and antipodal points', () => { + const acrossDateline = greatCircleKm( + { lat: 0, lon: 179 }, + { lat: 0, lon: -179 }, + ); + assert.ok(Math.abs(acrossDateline - 222.39) < 0.01); + assert.ok( + Math.abs( + greatCircleKm({ lat: 0, lon: 0 }, { lat: 0, lon: 180 }) - + Math.PI * EARTH_RADIUS_KM, + ) < 1e-9, + ); +}); + +test('canonical points require explicit latitude and longitude keys', () => { + assert.equal(isGeoPoint({ lat: 10, lon: 80 }), true); + assert.equal(isGeoPoint({ latitude: 10, longitude: 80 }), false); + assert.equal(isGeoPoint([10, 80]), false); + assert.equal( + greatCircleKm( + { latitude: 10, longitude: 80 }, + { latitude: 20, longitude: 70 }, + ), + Infinity, + ); +}); + +test('latitude-before-longitude regression rejects a swapped Austin point', () => { + assert.ok(Number.isFinite(greatCircleKm(AUSTIN, HOUSTON))); + assert.equal( + greatCircleKm( + { lat: AUSTIN.lon, lon: AUSTIN.lat }, + { lat: HOUSTON.lon, lon: HOUSTON.lat }, + ), + Infinity, + ); +}); + +test('invalid, missing, non-finite, and out-of-range points return Infinity', () => { + for (const point of [ + null, + {}, + { lat: 0 }, + { lat: NaN, lon: 0 }, + { lat: 0, lon: Infinity }, + { lat: 90.0001, lon: 0 }, + { lat: 0, lon: -180.0001 }, + ]) { + assert.equal(greatCircleKm(point, AUSTIN), Infinity); + assert.equal(greatCircleMeters(AUSTIN, point), Infinity); + } +}); diff --git a/src/hud.js b/src/hud.js index 783e1d566..8190f4b1b 100644 --- a/src/hud.js +++ b/src/hud.js @@ -15,6 +15,7 @@ import { applicationServices } from './services/application.js'; */ import * as Cesium from 'cesium'; +import { greatCircleKm } from './geoDistance.js'; import { forward as toMGRS } from 'mgrs'; import { CITY_POIS } from './cityPois.js'; import { composeLocalityTag } from './hudLocality.js'; @@ -683,13 +684,7 @@ export class IntelHUD { * @returns {number} Distance in kilometers. */ _haversineKm(lat1, lon1, lat2, lon2) { - const toRad = (deg) => Cesium.Math.toRadians(deg); - const dLat = toRad(lat2 - lat1); - const dLon = toRad(lon2 - lon1); - const a = - Math.sin(dLat / 2) ** 2 + - Math.cos(toRad(lat1)) * Math.cos(toRad(lat2)) * Math.sin(dLon / 2) ** 2; - return 6371 * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); + return greatCircleKm({ lat: lat1, lon: lon1 }, { lat: lat2, lon: lon2 }); } /** diff --git a/src/layers/bikeshare/model.js b/src/layers/bikeshare/model.js index d501d073c..96a543232 100644 --- a/src/layers/bikeshare/model.js +++ b/src/layers/bikeshare/model.js @@ -5,6 +5,7 @@ import { parseStationStatus, toNonNegativeInteger, } from '../../sources/gbfsStations.js'; +import { greatCircleKm } from '../../geoDistance.js'; export function createModel({ state: layerState, services, parts, source }) { /** @@ -28,15 +29,7 @@ export function createModel({ state: layerState, services, parts, source }) { */ function haversineKm(aLat, aLon, bLat, bLon) { - const toRad = (value) => (value * Math.PI) / 180; - const dLat = toRad(bLat - aLat); - const dLon = toRad(bLon - aLon); - const p1 = toRad(aLat); - const p2 = toRad(bLat); - const a = - Math.sin(dLat / 2) ** 2 + - Math.cos(p1) * Math.cos(p2) * Math.sin(dLon / 2) ** 2; - return 6371 * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); + return greatCircleKm({ lat: aLat, lon: aLon }, { lat: bLat, lon: bLon }); } /** diff --git a/src/layers/cctv/model.js b/src/layers/cctv/model.js index bdf350b4e..3354e0d43 100644 --- a/src/layers/cctv/model.js +++ b/src/layers/cctv/model.js @@ -1,4 +1,5 @@ import * as Cesium from 'cesium'; +import { greatCircleKm } from '../../geoDistance.js'; import { staticFrameRefreshMs } from '../../data/cctvLod.js'; import { frameFetchDue, cardFetchPolicy } from '../../data/cctvCards.js'; import { DEFAULT_CAMERA_CALIBRATION } from './policy.js'; @@ -309,12 +310,7 @@ export function createModel({ state: layerState, services, parts, source }) { */ function haversineKm(lat1, lon1, lat2, lon2) { - const dLat = toRad(lat2 - lat1); - const dLon = toRad(lon2 - lon1); - const a = - Math.sin(dLat / 2) ** 2 + - Math.cos(toRad(lat1)) * Math.cos(toRad(lat2)) * Math.sin(dLon / 2) ** 2; - return 6371 * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); + return greatCircleKm({ lat: lat1, lon: lon1 }, { lat: lat2, lon: lon2 }); } /** diff --git a/src/locations.js b/src/locations.js index 1c5e94f2f..d17186905 100644 --- a/src/locations.js +++ b/src/locations.js @@ -3,6 +3,7 @@ import { isUnavailableCapability } from './sources/capability.js'; import { createOverpassFeatureSource } from './sources/overpassFeatures.js'; import { applicationServices } from './services/application.js'; import * as Cesium from 'cesium'; +import { greatCircleKm } from './geoDistance.js'; import { viewportBias, placesNearViewRecovery, @@ -648,18 +649,6 @@ export const PLACE_ANCHOR_OFFSET_RATIO = 0.15; */ const PLACE_FALLBACK_HALF_SPAN_KM = 20; -/** Great-circle distance in km (small enough here that the spherical model is fine). */ -function greatCircleKm(lat1, lng1, lat2, lng2) { - const dLat = Cesium.Math.toRadians(lat2 - lat1); - const dLng = Cesium.Math.toRadians(lng2 - lng1); - const a = - Math.sin(dLat / 2) ** 2 + - Math.cos(Cesium.Math.toRadians(lat1)) * - Math.cos(Cesium.Math.toRadians(lat2)) * - Math.sin(dLng / 2) ** 2; - return 6371 * 2 * Math.asin(Math.min(1, Math.sqrt(a))); -} - /** * Sanity-gate a geocode viewport before framing it. Pure — exported for unit tests. * @@ -694,10 +683,8 @@ export function placeFramingViewport( return viewport; const offsetKm = greatCircleKm( - anchorLat, - anchorLng, - metrics.centerLat, - metrics.centerLng, + { lat: anchorLat, lon: anchorLng }, + { lat: metrics.centerLat, lon: metrics.centerLng }, ); if (offsetKm <= metrics.spanKm * PLACE_ANCHOR_OFFSET_RATIO) return viewport; diff --git a/src/voice/gevActions.js b/src/voice/gevActions.js index db089cecf..26974ce60 100644 --- a/src/voice/gevActions.js +++ b/src/voice/gevActions.js @@ -11,6 +11,7 @@ import { defaultGeospatial } from '../search/defaults.js'; import { msg } from '../i18n/index.js'; import * as Cesium from 'cesium'; import { CITY_POIS } from '../cityPois.js'; +import { greatCircleKm } from '../geoDistance.js'; import { findPoiByName, flyToGlobeView, @@ -3742,8 +3743,10 @@ function centerMatchesSubject(center, subjectPosition) { const subjectLat = Cesium.Math.toDegrees(carto.latitude); const subjectLon = Cesium.Math.toDegrees(carto.longitude); return ( - haversineKm(subjectLat, subjectLon, center.lat, center.lon) <= - SUBJECT_CENTER_TOLERANCE_KM + greatCircleKm( + { lat: subjectLat, lon: subjectLon }, + { lat: center.lat, lon: center.lon }, + ) <= SUBJECT_CENTER_TOLERANCE_KM ); } @@ -4033,7 +4036,7 @@ function nearbyKnownLandmarks(latitude, longitude, cameraHeightM) { for (const [cityId, city] of Object.entries(CITY_POIS)) { for (let poiIndex = 0; poiIndex < city.pois.length; poiIndex++) { const poi = city.pois[poiIndex]; - const distanceKm = haversineKm(latitude, longitude, poi.lat, poi.lon); + const distanceKm = greatCircleKm({ lat: latitude, lon: longitude }, poi); if (distanceKm > maxDistanceKm) continue; matches.push({ name: poi.name, @@ -4048,20 +4051,6 @@ function nearbyKnownLandmarks(latitude, longitude, cameraHeightM) { return matches.sort((a, b) => a.distanceKm - b.distanceKm).slice(0, 5); } -function haversineKm(lat1, lon1, lat2, lon2) { - const toRad = (value) => Cesium.Math.toRadians(value); - const radiusKm = 6371.0088; - const dLat = toRad(lat2 - lat1); - const dLon = toRad(lon2 - lon1); - const a = - Math.sin(dLat / 2) * Math.sin(dLat / 2) + - Math.cos(toRad(lat1)) * - Math.cos(toRad(lat2)) * - Math.sin(dLon / 2) * - Math.sin(dLon / 2); - return 2 * radiusKm * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); -} - function coarseBasemapPlace( viewScale, latitude,