impl/python.py

import math
from typing import Tuple

from .geo_bounding_box_types import BoundingBox

# 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).
#
# The trigonometry is written out below instead of calling math.sin and
# friends, because the platform maths library may differ in the last bit
# between Python, JavaScript and Rust. Only correctly-rounded IEEE 754
# operations (+ - * / sqrt floor) are used, 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

S3, S5, S7, S9, S11 = -1 / 6, 1 / 120, -1 / 5040, 1 / 362880, -1 / 39916800
S13, S15, S17 = 1 / 6227020800, -1 / 1307674368000, 1 / 355687428096000
C2, C4, C6, C8, C10 = -1 / 2, 1 / 24, -1 / 720, 1 / 40320, -1 / 3628800
C12, C14, C16, C18 = 1 / 479001600, -1 / 87178291200, 1 / 20922789888000, -1 / 6402373705728000
A3, A5, A7, A9, A11 = -1 / 3, 1 / 5, -1 / 7, 1 / 9, -1 / 11
A13, A15, A17, A19, A21 = 1 / 13, -1 / 15, 1 / 17, -1 / 19, 1 / 21


def _sin_small(r: float) -> float:
    s = r * r
    p = S17
    p = S15 + s * p
    p = S13 + s * p
    p = S11 + s * p
    p = S9 + s * p
    p = S7 + s * p
    p = S5 + s * p
    p = S3 + s * p
    return r + r * s * p


def _cos_small(r: float) -> float:
    s = r * r
    p = C18
    p = C16 + s * p
    p = C14 + s * p
    p = C12 + s * p
    p = C10 + s * p
    p = C8 + s * p
    p = C6 + s * p
    p = C4 + s * p
    p = C2 + s * p
    return 1 + s * p


def _sin_cos(x: float) -> Tuple[float, float]:
    k = math.floor(x / HALF_PI + 0.5)
    r = x - float(k) * HALF_PI
    sr = _sin_small(r)
    cr = _cos_small(r)
    quadrant = k % 4
    if quadrant == 0:
        return sr, cr
    if quadrant == 1:
        return cr, -sr
    if quadrant == 2:
        return -sr, -cr
    return -cr, sr


def _atan_unit(u: float) -> float:
    v = u
    v = v / (1 + math.sqrt(1 + v * v))
    v = v / (1 + math.sqrt(1 + v * v))
    v = v / (1 + math.sqrt(1 + v * v))
    s = v * v
    p = A21
    p = A19 + s * p
    p = A17 + s * p
    p = A15 + s * p
    p = A13 + s * p
    p = A11 + s * p
    p = A9 + s * p
    p = A7 + s * p
    p = A5 + s * p
    p = A3 + s * p
    return 8 * (v + v * s * p)


def _asin_positive(x: float) -> float:
    # (1 - x)(1 + x) keeps precision near 1.
    if x >= 1:
        return HALF_PI
    t = x / math.sqrt((1 - x) * (1 + x))
    return HALF_PI - _atan_unit(1 / t) if t > 1 else _atan_unit(t)


# 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)[0]
        cos_lat = _sin_cos(lat * DEGREES)[1]
        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),
    )