geo.bounding-box

Imports name this capability’s declared dependencies, which fune builds next to it in your project; each one links to its page.

use super::funejson::Value;  ← the fune runtime: the JSON value the test vectors use; fune build keeps it only where a signature takes one

// 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 f64::sin and
// friends, because the platform maths library may differ in the last bit
// between Rust, JavaScript and Python. 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)).
const EARTH_RADIUS_METRES: f64 = 6371008.8;
/// Ten millionths of a degree, about a centimetre.
const SCALE: f64 = 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.
const SNAP: f64 = 0.000001;

const PI: f64 = 3.141592653589793;
const HALF_PI: f64 = PI / 2.0;
const DEGREES: f64 = PI / 180.0;

const S3: f64 = -1.0 / 6.0;
const S5: f64 = 1.0 / 120.0;
const S7: f64 = -1.0 / 5040.0;
const S9: f64 = 1.0 / 362880.0;
const S11: f64 = -1.0 / 39916800.0;
const S13: f64 = 1.0 / 6227020800.0;
const S15: f64 = -1.0 / 1307674368000.0;
const S17: f64 = 1.0 / 355687428096000.0;
const C2: f64 = -1.0 / 2.0;
const C4: f64 = 1.0 / 24.0;
const C6: f64 = -1.0 / 720.0;
const C8: f64 = 1.0 / 40320.0;
const C10: f64 = -1.0 / 3628800.0;
const C12: f64 = 1.0 / 479001600.0;
const C14: f64 = -1.0 / 87178291200.0;
const C16: f64 = 1.0 / 20922789888000.0;
const C18: f64 = -1.0 / 6402373705728000.0;
const A3: f64 = -1.0 / 3.0;
const A5: f64 = 1.0 / 5.0;
const A7: f64 = -1.0 / 7.0;
const A9: f64 = 1.0 / 9.0;
const A11: f64 = -1.0 / 11.0;
const A13: f64 = 1.0 / 13.0;
const A15: f64 = -1.0 / 15.0;
const A17: f64 = 1.0 / 17.0;
const A19: f64 = -1.0 / 19.0;
const A21: f64 = 1.0 / 21.0;

fn sin_small(r: f64) -> f64 {
    let s = r * r;
    let mut 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;
    r + r * s * p
}

fn cos_small(r: f64) -> f64 {
    let s = r * r;
    let mut 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;
    1.0 + s * p
}

fn sin_cos(x: f64) -> (f64, f64) {
    let k = (x / HALF_PI + 0.5).floor();
    let r = x - k * HALF_PI;
    let sr = sin_small(r);
    let cr = cos_small(r);
    match (k as i64).rem_euclid(4) {
        0 => (sr, cr),
        1 => (cr, -sr),
        2 => (-sr, -cr),
        _ => (-cr, sr),
    }
}

fn atan_unit(u: f64) -> f64 {
    let mut v = u;
    v = v / (1.0 + (1.0 + v * v).sqrt());
    v = v / (1.0 + (1.0 + v * v).sqrt());
    v = v / (1.0 + (1.0 + v * v).sqrt());
    let s = v * v;
    let mut 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;
    8.0 * (v + v * s * p)
}

/// asin for 0 <= x, via atan; (1 - x)(1 + x) keeps precision near 1.
fn asin_positive(x: f64) -> f64 {
    if x >= 1.0 {
        return HALF_PI;
    }
    let t = x / ((1.0 - x) * (1.0 + x)).sqrt();
    if t > 1.0 {
        HALF_PI - atan_unit(1.0 / t)
    } 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.
fn round_down(x: f64) -> f64 {
    (x * SCALE + SNAP).floor() / SCALE + 0.0
}

fn round_up(x: f64) -> f64 {
    (x * SCALE - SNAP).ceil() / SCALE + 0.0
}

fn check_number(value: f64, what: &str) {
    if !value.is_finite() {
        panic!("{} must be a finite number", what);
    }
}

/// The box containing every point within `distance_metres` of (`lat`, `lng`).
///
/// # Panics
/// Panics if a coordinate is out of range or the distance is negative or not finite.
pub fn bounding_box(lat: f64, lng: f64, distance_metres: f64) -> BoundingBox {
    check_number(lat, "latitude");
    check_number(lng, "longitude");
    check_number(distance_metres, "distance");
    if lat < -90.0 || lat > 90.0 {
        panic!("latitude must be between -90 and 90 degrees");
    }
    if lng < -180.0 || lng > 180.0 {
        panic!("longitude must be between -180 and 180 degrees");
    }
    if distance_metres < 0.0 {
        panic!("distance must be 0 or more metres");
    }

    let r = distance_metres / EARTH_RADIUS_METRES;
    let r_degrees = r / DEGREES;
    let mut min_lat = lat - r_degrees;
    let mut max_lat = lat + r_degrees;
    let mut min_lng;
    let mut max_lng;

    if min_lat > -90.0 && max_lat < 90.0 {
        let sin_r = sin_cos(r).0;
        let cos_lat = sin_cos(lat * DEGREES).1;
        let 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.0 {
            min_lng += 360.0;
        }
        if max_lng > 180.0 {
            max_lng -= 360.0;
        }
    } else {
        // The circle reaches a pole, so it covers every longitude.
        if min_lat < -90.0 {
            min_lat = -90.0;
        }
        if max_lat > 90.0 {
            max_lat = 90.0;
        }
        min_lng = -180.0;
        max_lng = 180.0;
    }

    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),
    }
}

pub fn bounding_box_to_value(b: &BoundingBox) -> Value {
    Value::obj(vec![
        ("minLat", Value::Float(b.min_lat)),
        ("minLng", Value::Float(b.min_lng)),
        ("maxLat", Value::Float(b.max_lat)),
        ("maxLng", Value::Float(b.max_lng)),
    ])
}

pub fn fune_vector(args: &[Value]) -> Value {
    // Refuse what the typed signature cannot hold, with the wording TypeScript
    // and Python use, rather than let the conversion below quietly change it.
    if !matches!(args[0], Value::Int(_) | Value::Float(_)) {
        panic!("latitude must be a finite number");
    }
    if !matches!(args[1], Value::Int(_) | Value::Float(_)) {
        panic!("longitude must be a finite number");
    }
    if !matches!(args[2], Value::Int(_) | Value::Float(_)) {
        panic!("distance must be a finite number");
    }
    bounding_box_to_value(&bounding_box(
        args[0].as_f64(),
        args[1].as_f64(),
        args[2].as_f64(),
    ))
}