diff --git a/src/tools/area.js b/src/tools/area.js index e494c6782..41b417482 100644 --- a/src/tools/area.js +++ b/src/tools/area.js @@ -82,10 +82,18 @@ export async function resolveArea(area, { services, signal } = {}) { }; } +/** + * The bounding box of a spherical cap. When the cap reaches a pole it spans + * every longitude; otherwise the longitude span is the exact tangent bound, + * asin(sin(d) / cos(lat)) for angular radius d. + */ function circleArea(lat, lon, radiusKm) { - const latSpan = (radiusKm / EARTH_RADIUS_KM) * (180 / Math.PI); - const cosLat = Math.cos((lat * Math.PI) / 180); - const lonSpan = cosLat < 1e-6 ? 180 : Math.min(180, latSpan / cosLat); + const angular = radiusKm / EARTH_RADIUS_KM; + const latSpan = angular * (180 / Math.PI); + const reachesPole = lat + latSpan >= 90 || lat - latSpan <= -90; + const ratio = Math.sin(angular) / Math.cos((lat * Math.PI) / 180); + const lonSpan = + reachesPole || ratio >= 1 ? 180 : Math.asin(ratio) * (180 / Math.PI); const wrap = (value) => ((((value + 180) % 360) + 360) % 360) - 180; return { label: `${radiusKm} km around ${lat.toFixed(3)}, ${lon.toFixed(3)}`, diff --git a/src/tools/area.test.mjs b/src/tools/area.test.mjs index 99df85684..d22611282 100644 --- a/src/tools/area.test.mjs +++ b/src/tools/area.test.mjs @@ -117,3 +117,33 @@ test('a named place centers on its own point, not the middle of its bounds', asy const box = await resolveArea({ bbox: [0, 0, 10, 10] }); assert.deepEqual(areaCenter(box), { lat: 5, lon: 5 }); }); + +test('radius areas include every point within the radius, near the poles too', async () => { + const polar = await resolveArea({ lat: 80, lon: 0, radius_km: 1500 }); + const across = { lat: 89, lon: 180 }; + assert.ok(distanceKm({ lat: 80, lon: 0 }, across) < 1500); + assert.ok(areaContains(polar, across)); + assert.equal(polar.west, -180); + assert.equal(polar.east, 180); + assert.equal(polar.north, 90); + // Away from the poles, sample the circle's edge and check each point is kept. + const center = { lat: 60, lon: 20 }; + const wide = await resolveArea({ ...center, radius_km: 2000 }); + assert.ok(wide.east - wide.west < 360); + for (let bearing = 0; bearing < 360; bearing += 5) { + const d = 1999 / 6371.0088; + const b = (bearing * Math.PI) / 180; + const lat1 = (center.lat * Math.PI) / 180; + const lat2 = Math.asin( + Math.sin(lat1) * Math.cos(d) + Math.cos(lat1) * Math.sin(d) * Math.cos(b), + ); + const lon2 = + (center.lon * Math.PI) / 180 + + Math.atan2( + Math.sin(b) * Math.sin(d) * Math.cos(lat1), + Math.cos(d) - Math.sin(lat1) * Math.sin(lat2), + ); + const point = { lat: (lat2 * 180) / Math.PI, lon: (lon2 * 180) / Math.PI }; + assert.ok(areaContains(wide, point), `bearing ${bearing}`); + } +});