geo.bounding-box

export function boundingBox(lat: number, lng: number, distanceMetres: number): BoundingBox

/** Degrees, rounded outward to 7 decimals so the box always contains the circle. */
export interface BoundingBox {
  readonly minLat: number;
  readonly minLng: number;
  readonly maxLat: number;
  readonly maxLng: number;
}

import { boundingBox } from "#fune/geo.bounding-box@^2";

import { atan } from "./math_atan.ts";  ← from math.atan ^1.0.0 · built alongside by fune
import { sinCos } from "./math_sin_cos.ts";  ← from math.sin-cos ^1.0.0 · built alongside by fune
import { type BoundingBox } from "./geo_bounding_box_types.ts";

/**
 * The latitude/longitude box that contains every point within a distance of a
 * centre, on a sphere of the IUGG mean Earth radius, by the method of
 * J. P. Matuschek, "Finding Points Within a Distance of a Latitude/Longitude
 * Using Bounding Coordinates" (http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates).
 *
 * Sine, cosine and arctangent come from math.sin-cos and math.atan instead of
 * Math.sin and friends, because the platform maths library may differ in the
 * last bit between JavaScript, Python and Rust. The shared helpers use only
 * + - * /, and the rest here only + - * / sqrt floor ceil, all correctly
 * rounded, so every language computes the same doubles before rounding.
 */

// IUGG mean radius R1 of the GRS 80 ellipsoid (Moritz, Journal of Geodesy 74 (2000)).
const EARTH_RADIUS_METRES = 6371008.8;
// Ten millionths of a degree, about a centimetre.
const SCALE = 10000000;
// Values within a millionth of a grid step of a 7-decimal value are taken to
// be that value, so a centre of 51.5074 stays 51.5074 instead of becoming
// 51.5073999 because its double sits a hair below. That is about 1e-13
// degrees, far below anything the box is used for.
const SNAP = 0.000001;

const PI = 3.141592653589793;
const HALF_PI = PI / 2;
const DEGREES = PI / 180;

// asin for 0 <= x, on math.atan: asin x = atan(x / sqrt(1 - x^2)), with
// (1 - x)(1 + x) to keep precision near 1. There is no shared asin; x >= 1
// (reachable only by rounding, just short of a pole) is a quarter turn.
function asinPositive(x: number): number {
  if (x >= 1) return HALF_PI;
  return atan(x / Math.sqrt((1 - x) * (1 + x)));
}

// Outward rounding: a lower bound goes down and an upper bound goes up, so
// rounding can only grow the box, never cut the circle.
function roundDown(x: number): number {
  return Math.floor(x * SCALE + SNAP) / SCALE + 0;
}

function roundUp(x: number): number {
  return Math.ceil(x * SCALE - SNAP) / SCALE + 0;
}

function checkNumber(value: number, what: string): void {
  if (typeof value !== "number" || !Number.isFinite(value)) {
    throw new TypeError(`${what} must be a finite number`);
  }
}

export function boundingBox(lat: number, lng: number, distanceMetres: number): BoundingBox {
  checkNumber(lat, "latitude");
  checkNumber(lng, "longitude");
  checkNumber(distanceMetres, "distance");
  if (lat < -90 || lat > 90) throw new RangeError("latitude must be between -90 and 90 degrees");
  if (lng < -180 || lng > 180) throw new RangeError("longitude must be between -180 and 180 degrees");
  if (distanceMetres < 0) throw new RangeError("distance must be 0 or more metres");

  const r = distanceMetres / EARTH_RADIUS_METRES;
  const rDegrees = r / DEGREES;
  let minLat = lat - rDegrees;
  let maxLat = lat + rDegrees;
  let minLng: number;
  let maxLng: number;

  if (minLat > -90 && maxLat < 90) {
    const sinR = sinCos(r).sin;
    const cosLat = sinCos(lat * DEGREES).cos;
    const dLng = asinPositive(sinR / cosLat) / DEGREES;
    minLng = lng - dLng;
    maxLng = lng + dLng;
    // Past the antimeridian the bound comes round the other side, leaving
    // minLng > maxLng: the box is the two strips either side of 180.
    if (minLng < -180) minLng += 360;
    if (maxLng > 180) maxLng -= 360;
  } else {
    // The circle reaches a pole, so it covers every longitude.
    if (minLat < -90) minLat = -90;
    if (maxLat > 90) maxLat = 90;
    minLng = -180;
    maxLng = 180;
  }

  return { minLat: roundDown(minLat), minLng: roundDown(minLng), maxLat: roundUp(maxLat), maxLng: roundUp(maxLng) };
}

Install

fune build

With that line in your source, in a TypeScript project (language typescript in fune.project), fune build resolves it and its 2 dependencies, pins them in fune.lock, downloads only the TypeScript package of each, and builds the code above into your project’s .fune/build, one readable file per capability with a header linking back here. Or pin a range in fune.project and build in one step:

fune add geo.bounding-box

The manifest, vectors and README with only the TypeScript implementation. Install it without the registry with fune add ./geo.bounding-box-2.0.0-typescript.fune, or fetch it from a terminal with fune pull geo.bounding-box@2.0.0:typescript.

The whole function, every language, is one file too: geo.bounding-box-2.0.0.fune, 22,478 bytes, sha256 7e79132161c5e6092da15d9f79b040dfba351e8cb429691518d00652677f0457. It installs into a project of any language.

// fune: before geo.bounding-box

// fune: after geo.bounding-box

// fune: replace math.atan in geo.bounding-box
// fune: replace math.sin-cos in geo.bounding-box

// fune: step geo.bounding-box after <n|label>

More from the author

J. P. Matuschek, "Finding Points Within a Distance of a Latitude/Longitude Using Bounding Coordinates", http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates. The Earth is a sphere of radius 6,371,008.8 m (IUGG mean radius, as in `geo.distance`). The angular radius is r = d / R; latitudes are centre ± r; the longitude half-width is asin(sin r / cos lat), which is wider than r away from the equator (about twice as wide at 60°).

## Edge cases

- **Poles.** If the circle reaches or crosses a pole, it covers every longitude: the box is clamped to ±90 and runs from -180 to 180. Without this check the asin formula is asked for asin of more than 1 and returns NaN. - **The antimeridian.** A bound that passes ±180 comes round the other side, so minLng > maxLng. That is the signal that the box is two strips, minLng..180 and -180..maxLng; a query must OR them rather than use BETWEEN. - **Zero distance** is the point itself. **Negative distance** is an error. - A radius larger than the Earth gives the whole globe.

## Rounding outward

Bounds are rounded to 7 decimal places (about 1 cm), outward: the minimums down and the maximums up, so rounding can only grow the box and never cut off part of the circle. A value within a millionth of a step of a 7-decimal value (about 1e-13 degrees) is taken to be that value, so a centre of 51.5074 with distance 0 stays 51.5074 rather than becoming 51.5073999 because its binary value sits a hair below.

`math.round-float` is not used: it rounds to the nearest, and a box has to round each bound in its own direction, with the snap above. Those two one-line functions (floor or ceil of x × 10^7 ± 10^-6) stay here.

## What changed in 2.0.0

- **Sine, cosine and arctangent come from `math.sin-cos` and `math.atan`** instead of private copies. The registry has no asin, so a two-line asin for 0 ≤ x stays here, built on `math.atan`: asin x = atan(x / √((1 − x)(1 + x))), and π/2 when rounding pushes x to 1 or more just short of a pole. - **The unrounded bounds can differ from 1.x in the last bits** (for about one input in twelve), because the shared helpers reduce angles and evaluate atan differently. A difference that small changes a 7-decimal bound only when the bound lies within about 1e-13 degrees of a grid step, so no sampled box changed (none of 10,500), but it is not impossible, which is why this is a new major version. - **Vectors.** Every 1.x vector keeps its expected answer; each was re-derived from an independent 60-digit reference (own series for π, sin, cos and atan, asin through atan, decimal square root) with the outward rounding and snap applied to the true bounds, and matches. There is no rounding-change vector: the outward rounding is unchanged, and a trig-only change that crosses a grid step would need a bound within about 1e-13 degrees of it, roughly one chance in ten million per bound, too rare to construct from realistic coordinates. The vectors now compare exactly (`floats exact`): the three languages return the same doubles.

## Accuracy

Against the 60-digit reference, with the inputs taken as the exact doubles given, on 4,500 random centres within 80° of the equator and radii from 1 m to 20,000 km, the unrounded bounds are within 1.1e-13 degrees of the true bounds and the longitude half-width within 8 units in the last place. Within a fraction of a degree of a pole the problem is ill-conditioned (asin of nearly 1, the cosine of a latitude near 90°), and the error grows to about 3e-12 degrees, still some 30,000 times smaller than a grid step. 1.x had the same bounds. Every one of the 10,500 sampled boxes, rounded outward, equalled the true box rounded outward.

## Why the answers are identical in every language

The trigonometry comes from `math.sin-cos` and `math.atan`, which use only correctly rounded IEEE 754 operations (+ − × ÷) in the same order in all three languages, instead of the platform's Math.sin, which may differ in the last bit and could tip a value across a rounding step. The rest (sqrt, floor, ceil) is correctly rounded too, so TypeScript, Python and Rust compute the same doubles before and after rounding.

Files