impl/typescript.ts

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

import { atan } from "./math_atan.ts";  ← from math.atan ^1.0.0 · built alongside by fune
import { sinCos } from "./math_sin_cos.ts";  ← from math.sin-cos ^1.0.0 · built alongside by fune
import { type BoundingBox } from "./geo_bounding_box_types.ts";

/**
 * The latitude/longitude box that contains every point within a distance of a
 * centre, on a sphere of the IUGG mean Earth radius, 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 JavaScript, Python and Rust. The shared helpers use only
 * + - * /, and the rest here only + - * / sqrt floor ceil, all correctly
 * rounded, so every language computes the same doubles before rounding.
 */

// IUGG mean radius R1 of the GRS 80 ellipsoid (Moritz, Journal of Geodesy 74 (2000)).
const EARTH_RADIUS_METRES = 6371008.8;
// Ten millionths of a degree, about a centimetre.
const SCALE = 10000000;
// Values within a millionth of a grid step of a 7-decimal value are taken to
// be that value, so a centre of 51.5074 stays 51.5074 instead of becoming
// 51.5073999 because its double sits a hair below. That is about 1e-13
// degrees, far below anything the box is used for.
const SNAP = 0.000001;

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

// 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.
function asinPositive(x: number): number {
  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.
function roundDown(x: number): number {
  return Math.floor(x * SCALE + SNAP) / SCALE + 0;
}

function roundUp(x: number): number {
  return Math.ceil(x * SCALE - SNAP) / SCALE + 0;
}

function checkNumber(value: number, what: string): void {
  if (typeof value !== "number" || !Number.isFinite(value)) {
    throw new TypeError(`${what} must be a finite number`);
  }
}

export function boundingBox(lat: number, lng: number, distanceMetres: number): BoundingBox {
  checkNumber(lat, "latitude");
  checkNumber(lng, "longitude");
  checkNumber(distanceMetres, "distance");
  if (lat < -90 || lat > 90) throw new RangeError("latitude must be between -90 and 90 degrees");
  if (lng < -180 || lng > 180) throw new RangeError("longitude must be between -180 and 180 degrees");
  if (distanceMetres < 0) throw new RangeError("distance must be 0 or more metres");

  const r = distanceMetres / EARTH_RADIUS_METRES;
  const rDegrees = r / DEGREES;
  let minLat = lat - rDegrees;
  let maxLat = lat + rDegrees;
  let minLng: number;
  let maxLng: number;

  if (minLat > -90 && maxLat < 90) {
    const sinR = sinCos(r).sin;
    const cosLat = sinCos(lat * DEGREES).cos;
    const dLng = asinPositive(sinR / cosLat) / DEGREES;
    minLng = lng - dLng;
    maxLng = lng + dLng;
    // Past the antimeridian the bound comes round the other side, leaving
    // minLng > maxLng: the box is the two strips either side of 180.
    if (minLng < -180) minLng += 360;
    if (maxLng > 180) maxLng -= 360;
  } else {
    // The circle reaches a pole, so it covers every longitude.
    if (minLat < -90) minLat = -90;
    if (maxLat > 90) maxLat = 90;
    minLng = -180;
    maxLng = 180;
  }

  return { minLat: roundDown(minLat), minLng: roundDown(minLng), maxLat: roundUp(maxLat), maxLng: roundUp(maxLng) };
}