geo_bounding_box.rb

class GeoBoundingBox

  attr_accessor :lat, :lon, :distance

  WGS84_MAJOR_AXIS = 6378137
  WGS84_MINOR_AXIS = 6356752.3142
  WGS84_FLATTENING = 1 / 298.257223563

  def initialize(lat, lon, distance)
    @lat = lat
    @lon = lon
    @distance = distance
  end

  def north_east
    point_from_distance_with_bearing(@lat, @lon, @distance / 2, 45)
  end

  def south_west
    point_from_distance_with_bearing(@lat, @lon, @distance / 2, 225)
  end

  def to_radians(degrees)
    degrees * Math::PI / 180
  end

  def to_degrees(radians)
    radians * 180 / Math::PI
  end

  private

    # implementation of Vincenty formula, from:
    # http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html
    # distance is in m
    # bearing is in decimal degrees
    def point_from_distance_with_bearing(lat, lon, distance, bearing)

      a = WGS84_MAJOR_AXIS
      b = WGS84_MINOR_AXIS
      f = WGS84_FLATTENING

      s = distance;
      alpha1 = to_radians(bearing)
      sinAlpha1 = Math.sin(alpha1);
      cosAlpha1 = Math.cos(alpha1);

      tanU1 = (1 - f) * Math.tan(to_radians(lat))
      cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1))
      sinU1 = tanU1 * cosU1
      sigma1 = Math.atan2(tanU1, cosAlpha1)
      sinAlpha = cosU1 * sinAlpha1
      cosSqAlpha = 1 - sinAlpha * sinAlpha
      uSq = cosSqAlpha * (a * a - b * b) / (b * b)
      a2 = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
      b2 = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))

      sigma = s / (b * a2)
      sigmaP = 2 * Math::PI

      while ((sigma - sigmaP).abs > 1e-12) do
        cos2SigmaM = Math.cos(2 * sigma1 + sigma)
        sinSigma = Math.sin(sigma)
        cosSigma = Math.cos(sigma)
        deltaSigma = b2 * sinSigma * (cos2SigmaM + b2 / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)\
          - b2 / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)))
        sigmaP = sigma
        sigma = s / (b * a2) + deltaSigma
      end

      tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1
      lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,\
        (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp))

      lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1)
      c = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
      l = lambda - (1 - c) * f * sinAlpha * (sigma + c * sinSigma * (cos2SigmaM + c * cosSigma*(-1 + 2 * cos2SigmaM * cos2SigmaM)))
      lon2 = (to_radians(lon) + l + 3 * Math::PI) % (2 * Math::PI) - Math::PI # normalise to -180...+180

      revAz = Math.atan2(sinAlpha, -tmp) # final bearing, if required

      [ to_degrees(lat2), to_degrees(lon2) ]
    end

end