Geo Bunding Box

geoBoundingBox.ts

/**
 * Utility functions:
 * - getBoundingBox: Gets a box surrounding all points that are a certain distance from provided point on the Earth sphere.
 *                   Take care when longitude crosses the 180th meridian. In such case, the minLongitude (e.g., 175) 
 *                   will be greater than the maxLongitude (e.g. -175).
 *                   Credit to Jan Philip Matuschek for the mathematical formulas - http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
 * - getDistanceInMeters: Gets a distance between 2 points on the Earth sphere.
 */
export type BoundingBox = { minLatitude: number; maxLatitude: number; minLongitude: number; maxLongitude: number }

const degToRad = (deg: number): number => (deg * Math.PI) / 180.0
const radToDeg = (rad: number): number => (rad * 180.0) / Math.PI

const EARTH_RADIUS_METERS =  6_356_752.3

const MIN_LATITUDE_RAD = -Math.PI / 2
const MAX_LATITUDE_RAD = Math.PI / 2
const MIN_LONGITUDE_RAD = -Math.PI
const MAX_LONGITUDE_RAD = Math.PI

export const getDistanceInMeters = (
  latitude1: number,
  longitude1: number,
  latitude2: number,
  longitude2: number
): number => {
  const latitudeRad1 = degToRad(latitude1)
  const longitudeRad1 = degToRad(longitude1)
  const latitudeRad2 = degToRad(latitude2)
  const longitudeRad2 = degToRad(longitude2)

  return (
    Math.acos(
      Math.sin(latitudeRad1) * Math.sin(latitudeRad2) +
        Math.cos(latitudeRad1) * Math.cos(latitudeRad2) * Math.cos(longitudeRad1 - longitudeRad2)
    ) * EARTH_RADIUS_METERS
  )
}

export const getBoundingBox = (latitude: number, longitude: number, distanceInMeters: number): BoundingBox => {
  const latitudeRad = degToRad(latitude)
  const longitudeRad = degToRad(longitude)

  // angular distance in radians on a great circle
  const distanceRad = distanceInMeters / EARTH_RADIUS_METERS

  let minLatitudeRad = latitudeRad - distanceRad
  let maxLatitudeRad = latitudeRad + distanceRad

  let minLongitudeRad: number
  let maxLongitudeRad: number
  if (minLatitudeRad > MIN_LATITUDE_RAD && maxLatitudeRad < MAX_LATITUDE_RAD) {
    const deltaLongitudeRad = Math.asin(Math.sin(distanceRad) / Math.cos(latitudeRad))

    minLongitudeRad = longitudeRad - deltaLongitudeRad
    if (minLongitudeRad < MIN_LONGITUDE_RAD) {
      // Overflow over the 180th meridian "to the other end of the map"
      minLongitudeRad += 2.0 * Math.PI
    }

    maxLongitudeRad = longitudeRad + deltaLongitudeRad
    if (maxLongitudeRad > MAX_LONGITUDE_RAD) {
      // Overflow over the 180th meridian "to the other end of the map"
      maxLongitudeRad -= 2.0 * Math.PI
    }
  } else {
    // Stopping at north or south pole
    minLatitudeRad = Math.max(minLatitudeRad, MIN_LATITUDE_RAD)
    maxLatitudeRad = Math.min(maxLatitudeRad, MAX_LATITUDE_RAD)
    // If earth pole is in the bounding box, we take the whole horizontal slice.
    minLongitudeRad = MIN_LONGITUDE_RAD
    maxLongitudeRad = MAX_LONGITUDE_RAD
  }

  return {
    minLatitude: radToDeg(minLatitudeRad),
    maxLatitude: radToDeg(maxLatitudeRad),
    minLongitude: radToDeg(minLongitudeRad),
    maxLongitude: radToDeg(maxLongitudeRad),
  }
}