Functional Weave
Code in Python

geo.distance@2.0.0

impl/python.py

2,856 bytes · the Python implementation · view raw

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 .math_atan import atan  ← from math.atan ^1.0.0 · built alongside by fune
from .math_round_float import round_float  ← from math.round-float ^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

# Sine, cosine and arctangent come from math.sin-cos and math.atan, not
# math.sin or math.atan2: those come from the platform maths library, which may
# differ in the last bit between Python, JavaScript and Rust, and a value next
# to a rounding boundary can then round differently. With the shared helpers
# the unrounded distance is the same double everywhere, and math.round-float
# rounds it the same way everywhere.

# IUGG mean radius R1 = (2a + b) / 3 of the GRS 80 ellipsoid (Moritz,
# "Geodetic Reference System 1980", Journal of Geodesy 74 (2000) 128-133).
EARTH_RADIUS_METRES = 6371008.8

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


def _atan2_positive(y: float, x: float) -> float:
    # atan2 for y >= 0 and x >= 0. x is 0 only for exactly antipodal points,
    # where y / x would divide by zero.
    if x == 0:
        return 0.0 if y == 0 else HALF_PI
    return atan(y / x)


def _is_number(value: object) -> bool:
    return isinstance(value, (int, float)) and not isinstance(value, bool) and math.isfinite(value)


def _check_latitude(value: float) -> None:
    if not _is_number(value):
        raise TypeError("latitude must be a finite number of degrees")
    if value < -90 or value > 90:
        raise ValueError("latitude must be between -90 and 90 degrees")


def _check_longitude(value: float) -> None:
    if not _is_number(value):
        raise TypeError("longitude must be a finite number of degrees")
    if value < -180 or value > 180:
        raise ValueError("longitude must be between -180 and 180 degrees")


def distance(from_lat: float, from_lng: float, to_lat: float, to_lng: float) -> float:
    """Great-circle distance in metres (haversine, IUGG mean radius), to the millimetre."""
    _check_latitude(from_lat)
    _check_longitude(from_lng)
    _check_latitude(to_lat)
    _check_longitude(to_lng)
    from_lat, from_lng, to_lat, to_lng = float(from_lat), float(from_lng), float(to_lat), float(to_lng)

    # A longitude difference of 359 degrees is a 1 degree step across the
    # antimeridian; sin^2 of the half-angle is the same either way.
    sin_half_lat = sin_cos(((to_lat - from_lat) * DEGREES) / 2).sin
    sin_half_lng = sin_cos(((to_lng - from_lng) * DEGREES) / 2).sin
    cos1 = sin_cos(from_lat * DEGREES).cos
    cos2 = sin_cos(to_lat * DEGREES).cos

    a = sin_half_lat * sin_half_lat + cos1 * cos2 * sin_half_lng * sin_half_lng
    if a < 0:
        a = 0.0
    if a > 1:
        a = 1.0
    # atan2 rather than asin(sqrt(a)): stays accurate for antipodal points too.
    c = 2 * _atan2_positive(math.sqrt(a), math.sqrt(1 - a))
    # Half away from zero on the exact value of the double (math.round-float).
    return round_float(EARTH_RADIUS_METRES * c, 3)