Imports name this capability’s declared dependencies, which fune builds next to it in your project; each one links to its page.
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),
)