geo.bounding-box

def bounding_box(lat: float, lng: float, distance_metres: float) -> BoundingBox

@dataclass(frozen=True)
class BoundingBox:
    """Degrees, rounded outward to 7 decimals so the box always contains the circle."""

    min_lat: float
    min_lng: float
    max_lat: float
    max_lng: float

from fune.geo.bounding_box import bounding_box  # geo.bounding-box@^2

import math

from .geo_bounding_box_types import BoundingBox
from .math_atan import atan  ← from math.atan ^1.0.0 · built alongside by fune
from .math_sin_cos import sin_cos  ← from math.sin-cos ^1.0.0 · built alongside by fune

# The latitude/longitude box that contains every point within a distance of a
# centre, 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 Python, JavaScript and Rust. The shared helpers use only
# + - * /, and the rest here only + - * / sqrt floor ceil, all correctly
# rounded, in the same order in all three.

# IUGG mean radius R1 of the GRS 80 ellipsoid (Moritz, Journal of Geodesy 74 (2000)).
EARTH_RADIUS_METRES = 6371008.8
# Ten millionths of a degree, about a centimetre.
SCALE = 10000000.0
# Within a millionth of a grid step of a 7-decimal value counts as that value,
# so a centre of 51.5074 stays 51.5074 rather than becoming 51.5073999.
SNAP = 0.000001

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


def _asin_positive(x: float) -> float:
    # 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.
    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.
def _round_down(x: float) -> float:
    return float(math.floor(x * SCALE + SNAP)) / SCALE + 0.0


def _round_up(x: float) -> float:
    return float(math.ceil(x * SCALE - SNAP)) / SCALE + 0.0


def _check_number(value: object, what: str) -> None:
    if isinstance(value, bool) or not isinstance(value, (int, float)) or not math.isfinite(value):
        raise TypeError("%s must be a finite number" % what)


def bounding_box(lat: float, lng: float, distance_metres: float) -> BoundingBox:
    """The box containing every point within distance_metres of (lat, lng)."""
    _check_number(lat, "latitude")
    _check_number(lng, "longitude")
    _check_number(distance_metres, "distance")
    if lat < -90 or lat > 90:
        raise ValueError("latitude must be between -90 and 90 degrees")
    if lng < -180 or lng > 180:
        raise ValueError("longitude must be between -180 and 180 degrees")
    if distance_metres < 0:
        raise ValueError("distance must be 0 or more metres")
    lat, lng, distance_metres = float(lat), float(lng), float(distance_metres)

    r = distance_metres / EARTH_RADIUS_METRES
    r_degrees = r / DEGREES
    min_lat = lat - r_degrees
    max_lat = lat + r_degrees

    if min_lat > -90 and max_lat < 90:
        sin_r = sin_cos(r).sin
        cos_lat = sin_cos(lat * DEGREES).cos
        d_lng = _asin_positive(sin_r / cos_lat) / DEGREES
        min_lng = lng - d_lng
        max_lng = lng + d_lng
        # Past the antimeridian the bound comes round the other side, leaving
        # min_lng > max_lng: the box is the two strips either side of 180.
        if min_lng < -180:
            min_lng += 360
        if max_lng > 180:
            max_lng -= 360
    else:
        # The circle reaches a pole, so it covers every longitude.
        if min_lat < -90:
            min_lat = -90.0
        if max_lat > 90:
            max_lat = 90.0
        min_lng = -180.0
        max_lng = 180.0

    return BoundingBox(
        min_lat=_round_down(min_lat),
        min_lng=_round_down(min_lng),
        max_lat=_round_up(max_lat),
        max_lng=_round_up(max_lng),
    )

Install

fune build

With that line in your source, in a Python project (language python in fune.project), fune build resolves it and its 2 dependencies, pins them in fune.lock, downloads only the Python 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 Python implementation. Install it without the registry with fune add ./geo.bounding-box-2.0.0-python.fune, or fetch it from a terminal with fune pull geo.bounding-box@2.0.0:python.

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