Functional Weave
Code in Python

geo.bounding-box@2.0.0

README.md

4,494 bytes · view raw

# geo.bounding-box

The smallest latitude/longitude box that contains every point within a given
distance of a centre. It is the cheap first filter for "find everything within
10 km": query the box with an index, then check the real distance
(`geo.distance`) on what comes back.

## Method

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.