Calculate distance between two points on a globe

From CodeCodex

Contents

Haversine Formula

The haversine formula is an equation important in navigation, giving great-circle distances between two points on a sphere from their longitudes and latitudes. It is a special case of a more general formula in spherical trigonometry, the law of haversines, relating the sides and angles of spherical "triangles".

Implementations

APEX Force.com


public class LocatorUtil {

private static Double EARTH_RADIUS = 6371.00; // Radius in Kilometers default
		
public static Double calculateDistance(Double lat1, Double lon1, Double lat2,   Double lon2){
  Double Radius = LocatorUtil.EARTH_RADIUS; //6371.00;
  Double dLat = LocatorUtil.toRadians(lat2-lat1);
  Double dLon = LocatorUtil.toRadians(lon2-lon1);            
  Double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
  Math.cos(LocatorUtil.toRadians(lat1)) *   Math.cos(LocatorUtil.toRadians(lat2)) *
  Math.sin(dLon/2) * Math.sin(dLon/2);
  Double c = 2 * Math.asin(Math.sqrt(a));
  return Radius * c;		
}
	
public static Double toRadians(Double degree){
   // Value degree * Pi/180
   Double res = degree * 3.1415926 / 180;
   return res;	
}
}

Java

import com.google.android.maps.GeoPoint;

public class DistanceCalculator {

   private double Radius;

   // R = earth's radius (mean radius = 6,371km)
   // Constructor
   DistanceCalculator(double R) {
      Radius = R;
   }

   public double CalculationByDistance(GeoPoint StartP, GeoPoint EndP) {
      double lat1 = StartP.getLatitudeE6()/1E6;
      double lat2 = EndP.getLatitudeE6()/1E6;
      double lon1 = StartP.getLongitudeE6()/1E6;
      double lon2 = EndP.getLongitudeE6()/1E6;
      double dLat = Math.toRadians(lat2-lat1);
      double dLon = Math.toRadians(lon2-lon1);
      double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
         Math.cos(Math.toRadians(lat1)) * Math.cos(Math.toRadians(lat2)) *
         Math.sin(dLon/2) * Math.sin(dLon/2);
      double c = 2 * Math.asin(Math.sqrt(a));
      return Radius * c;
   }
}


Objective C

// adapted from C++ code on this page
+ (CGFloat)directMetersFromCoordinate:(CLLocationCoordinate2D)from toCoordinate:(CLLocationCoordinate2D)to {
    
    static const double DEG_TO_RAD = 0.017453292519943295769236907684886;  
    static const double EARTH_RADIUS_IN_METERS = 6372797.560856;  
    
    double latitudeArc  = (from.latitude - to.latitude) * DEG_TO_RAD;
    double longitudeArc = (from.longitude - to.longitude) * DEG_TO_RAD;
    double latitudeH = sin(latitudeArc * 0.5);
    latitudeH *= latitudeH;
    double lontitudeH = sin(longitudeArc * 0.5);
    lontitudeH *= lontitudeH;
    double tmp = cos(from.latitude*DEG_TO_RAD) * cos(to.latitude*DEG_TO_RAD);
    return EARTH_RADIUS_IN_METERS * 2.0 * asin(sqrt(latitudeH + tmp*lontitudeH));
}

Perl

# $lat1 and $lon1 are the coordinates of the first point in radians
# $lat2 and $lon2 are the coordinates of the second point in radians
my $a = sin(($lat2 - $lat1)/2.0);
my $b = sin(($lon2 - $lon1)/2.0);
my $h = ($a*$a) + cos($lat1) * cos($lat2) * ($b*$b);
my $theta = 2 * asin(sqrt($h)); # distance in radians
# in order to find the distance, multiply $theta by the radius of the earth, e.g.
# $theta * 6,372.7976 = distance in kilometres (value from http://en.wikipedia.org/wiki/Earth_radius)

PostgreSQL PL/PERLU

create or replace function distance_in_km (numeric,numeric,numeric,numeric) returns numeric as
$body$
use strict;
use Math::Trig qw(great_circle_distance deg2rad);

my $lat1 = shift(@_);
my $lon1 = shift(@_);
my $lat2 = shift(@_);
my $lon2 = shift(@_);

# Notice the 90 - latitude: phi zero is at the North Pole.
sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
my @L = NESW( $lat1, $lon1 );
my @T = NESW( $lat2, $lon2 );
my $km = great_circle_distance(@L, @T, 6378);

return $km;
$body$
language plperlu strict security definer;

Query that returns nearest 10 locations within 25 km. Assumes the location is stored as two numeric fields, lat and lng.

SELECT *, (
	6371 * acos(cos(radians(search_lat)) * cos(radians(lat) ) * 
	cos(radians(lng) - radians(search_lng)) + sin(radians(search_lat)) * sin(radians(lat)))
) AS distance
FROM table
WHERE lat != search_lat AND lng != search_lng AND distance < 25
ORDER BY distance
FETCH 10 ONLY

JavaScript

var R = 6371; // km
var dLat = (lat2-lat1)*Math.PI/180;
var dLon = (lon2-lon1)*Math.PI/180; 
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
        Math.cos(lat1*Math.PI/180) * Math.cos(lat2*Math.PI/180) * 
        Math.sin(dLon/2) * Math.sin(dLon/2); 
var c = 2 * Math.asin(Math.sqrt(a)); 
var d = R * c;

C++

/// @brief The usual PI/180 constant
static const double DEG_TO_RAD = 0.017453292519943295769236907684886;
/// @brief Earth's quatratic mean radius for WGS-84
static const double EARTH_RADIUS_IN_METERS = 6372797.560856;

/** @brief Computes the arc, in radian, between two WGS-84 positions.
  *
  * The result is equal to <code>Distance(from,to)/EARTH_RADIUS_IN_METERS</code>
  *    <code>= 2*asin(sqrt(h(d/EARTH_RADIUS_IN_METERS )))</code>
  *
  * where:<ul>
  *    <li>d is the distance in meters between 'from' and 'to' positions.</li>
  *    <li>h is the haversine function: <code>h(x)=sin²(x/2)</code></li>
  * </ul>
  *
  * The haversine formula gives:
  *    <code>h(d/R) = h(from.lat-to.lat)+h(from.lon-to.lon)+cos(from.lat)*cos(to.lat)</code>
  *
  * @sa http://en.wikipedia.org/wiki/Law_of_haversines
  */
double ArcInRadians(const Position& from, const Position& to) {
    double latitudeArc  = (from.lat - to.lat) * DEG_TO_RAD;
    double longitudeArc = (from.lon - to.lon) * DEG_TO_RAD;
    double latitudeH = sin(latitudeArc * 0.5);
    latitudeH *= latitudeH;
    double lontitudeH = sin(longitudeArc * 0.5);
    lontitudeH *= lontitudeH;
    double tmp = cos(from.lat*DEG_TO_RAD) * cos(to.lat*DEG_TO_RAD);
    return 2.0 * asin(sqrt(latitudeH + tmp*lontitudeH));
}

/** @brief Computes the distance, in meters, between two WGS-84 positions.
  *
  * The result is equal to <code>EARTH_RADIUS_IN_METERS*ArcInRadians(from,to)</code>
  *
  * @sa ArcInRadians
  */
double DistanceInMeters(const Position& from, const Position& to) {
    return EARTH_RADIUS_IN_METERS*ArcInRadians(from, to);
}

Ruby

# haversine.rb
#
# haversine formula to compute the great circle distance between two points given their latitude and longitudes
#
# Copyright (C) 2008, 360VL, Inc
# Copyright (C) 2008, Landon Cox
#
# http://www.esawdust.com (Landon Cox)
# contact:
# http://www.esawdust.com/blog/businesscard/businesscard.html
#
# LICENSE: GNU Affero GPL v3
# The ruby implementation of the Haversine formula is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License version 3 as published by the Free Software Foundation. 
#
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public
# License version 3 for more details.  http://www.gnu.org/licenses/
#
# Landon Cox - 9/25/08
#
# Notes:
#
# translated into Ruby based on information contained in:
#   http://mathforum.org/library/drmath/view/51879.html  Doctors Rick and Peterson - 4/20/99
#   http://www.movable-type.co.uk/scripts/latlong.html
#   http://en.wikipedia.org/wiki/Haversine_formula
#
# This formula can compute accurate distances between two points given latitude and longitude, even for
# short distances.
 
# PI = 3.1415926535
RAD_PER_DEG = 0.017453293  #  PI/180

# the great circle distance d will be in whatever units R is in

Rmiles = 3956           # radius of the great circle in miles
Rkm = 6371              # radius in kilometers...some algorithms use 6367
Rfeet = Rmiles * 5282   # radius in feet
Rmeters = Rkm * 1000    # radius in meters

@distances = Hash.new   # this is global because if computing lots of track point distances, it didn't make
                        # sense to new a Hash each time over potentially 100's of thousands of points

=begin rdoc
  given two lat/lon points, compute the distance between the two points using the haversine formula
  the result will be a Hash of distances which are key'd by 'mi','km','ft', and 'm'
=end

def haversine_distance( lat1, lon1, lat2, lon2 )

  dlon = lon2 - lon1
  dlat = lat2 - lat1
 
  dlon_rad = dlon * RAD_PER_DEG
  dlat_rad = dlat * RAD_PER_DEG
 
  lat1_rad = lat1 * RAD_PER_DEG
  lon1_rad = lon1 * RAD_PER_DEG
 
  lat2_rad = lat2 * RAD_PER_DEG
  lon2_rad = lon2 * RAD_PER_DEG
 
  # puts "dlon: #{dlon}, dlon_rad: #{dlon_rad}, dlat: #{dlat}, dlat_rad: #{dlat_rad}"
 
  a = Math.sin(dlat_rad/2)**2 + Math.cos(lat1_rad) * Math.cos(lat2_rad) * Math.sin(dlon_rad/2)**2
  c = 2 * Math.asin( Math.sqrt(a))
 
  dMi = Rmiles * c          # delta between the two points in miles
  dKm = Rkm * c             # delta in kilometers
  dFeet = Rfeet * c         # delta in feet
  dMeters = Rmeters * c     # delta in meters
 
  @distances["mi"] = dMi
  @distances["km"] = dKm
  @distances["ft"] = dFeet
  @distances["m"] = dMeters
end

def test_haversine
 
  lon1 = -104.88544
  lat1 = 39.06546
 
  lon2 = -104.80
  lat2 = lat1
 
  haversine_distance( lat1, lon1, lat2, lon2 )
 
  puts "the distance from  #{lat1}, #{lon1} to #{lat2}, #{lon2} is:"
  puts "#{@distances['mi']} mi"
  puts "#{@distances['km']} km"
  puts "#{@distances['ft']} ft"
  puts "#{@distances['m']} m"
 
  if @distances['km'].to_s =~ /7\.376*/
    puts "Test: Success"
  else
    puts "Test: Failed"
  end
 
end

test_haversine

Python

#coding:UTF-8
"""
  Python implementation of Haversine formula
  Copyright (C) <2009>  Bartek Górny <bartek@gorny.edu.pl>

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

import math

def recalculate_coordinate(val,  _as=None):
  """
    Accepts a coordinate as a tuple (degree, minutes, seconds)
    You can give only one of them (e.g. only minutes as a floating point number) and it will be duly
    recalculated into degrees, minutes and seconds.
    Return value can be specified as 'deg', 'min' or 'sec'; default return value is a proper coordinate tuple.
  """
  deg,  min,  sec = val
  # pass outstanding values from right to left
  min = (min or 0) + int(sec) / 60
  sec = sec % 60
  deg = (deg or 0) + int(min) / 60
  min = min % 60
  # pass decimal part from left to right
  dfrac,  dint = math.modf(deg)
  min = min + dfrac * 60
  deg = dint
  mfrac,  mint = math.modf(min)
  sec = sec + mfrac * 60
  min = mint
  if _as:
    sec = sec + min * 60 + deg * 3600
    if _as == 'sec': return sec
    if _as == 'min': return sec / 60
    if _as == 'deg': return sec / 3600
  return deg,  min,  sec
      

def points2distance(start,  end):
  """
    Calculate distance (in kilometers) between two points given as (long, latt) pairs
    based on Haversine formula (http://en.wikipedia.org/wiki/Haversine_formula).
    Implementation inspired by JavaScript implementation from http://www.movable-type.co.uk/scripts/latlong.html
    Accepts coordinates as tuples (deg, min, sec), but coordinates can be given in any form - e.g.
    can specify only minutes:
    (0, 3133.9333, 0) 
    is interpreted as 
    (52.0, 13.0, 55.998000000008687)
    which, not accidentally, is the lattitude of Warsaw, Poland.
  """
  start_long = math.radians(recalculate_coordinate(start[0],  'deg'))
  start_latt = math.radians(recalculate_coordinate(start[1],  'deg'))
  end_long = math.radians(recalculate_coordinate(end[0],  'deg'))
  end_latt = math.radians(recalculate_coordinate(end[1],  'deg'))
  d_latt = end_latt - start_latt
  d_long = end_long - start_long
  a = math.sin(d_latt/2)**2 + math.cos(start_latt) * math.cos(end_latt) * math.sin(d_long/2)**2
  c = 2 * math.asin(math.sqrt(a))
  return 6371 * c


if __name__ == '__main__':
 warsaw = ((21,  0,  30),  (52, 13, 56))
 cracow = ((19, 56, 18),  (50, 3, 41))
 print points2distance(warsaw,  cracow)

MySQL

DELIMITER //
CREATE DEFINER=`someone`@`somehost` FUNCTION `GetDistance`(
 lat1  numeric (9,6),
 lon1  numeric (9,6),
 lat2  numeric (9,6),
 lon2  numeric (9,6)
) RETURNS decimal(10,5)
    READS SQL DATA
BEGIN
  DECLARE  x  decimal (20,10);
  DECLARE  pi  decimal (21,20);
  SET  pi = 3.14159265358979323846;
  SET  x = sin( lat1 * pi/180 ) * sin( lat2 * pi/180  ) + cos(
 lat1 *pi/180 ) * cos( lat2 * pi/180 ) * cos(  abs( (lon2 * pi/180) -
 (lon1 *pi/180) ) );
  SET  x = acos( x );
  RETURN  ( 1.852 * 60.0 * ((x/pi)*180) ) / 1.609344;
END //
Calculate the distance between 2 points with latitude & longitude in MySQL:

DISTANCE in Kilometres:

6371 * acos( cos( radians(LATITUDE1) ) * cos( radians( LATITUDE2 ) ) * cos( radians( LONGITUDE1 ) - radians(LONGITUDE2) ) + sin( radians(LATITUDE1) ) * sin( radians( LATITUDE2 ) ) ) AS DISTANCE  

DISTANCE in Miles:

3956 * acos( cos( radians(LATITUDE1) ) * cos( radians( LATITUDE2 ) ) * cos( radians( LONGITUDE1 ) - radians(LONGITUDE2) ) + sin( radians(LATITUDE1) ) * sin( radians( LATITUDE2 ) ) ) AS DISTANCE  


List all points in table having distance between a designated point (we use a random point - lat:45.20327, long:23.7806) less than 50 KM, with latitude & longitude, in MySQL (the table fields are coord_lat and coord_long):

List all having DISTANCE<50, in Kilometres (considered Earth radius 6371 KM):

SELECT denumire, (6371 * acos( cos( radians(45.20327) ) * cos( radians( coord_lat ) ) * cos( radians( 23.7806 ) - radians(coord_long) ) + sin( radians(45.20327) ) * sin( radians(coord_lat) ) )) AS distanta FROM obiective WHERE coord_lat<>'' AND coord_long<>'' having distanta<50 ORDER BY distanta desc

The above example was tested in MySQL 5.0.95 and 5.5.16 (Linux).

MSSQL

CREATE FUNCTION [dbo].[GetDistance]
(
      @lat1 Float(8),
      @long1 Float(8),
      @lat2 Float(8),
      @long2 Float(8)
)
RETURNS Float(8)
AS
BEGIN
      DECLARE @R Float(8);
      DECLARE @dLat Float(8);
      DECLARE @dLon Float(8);
      DECLARE @a Float(8);
      DECLARE @c Float(8);
      DECLARE @d Float(8);
 
      SET @R = 3960;
      SET @dLat = RADIANS(@lat2 - @lat1);
      SET @dLon = RADIANS(@long2 - @long1);
 
      SET @a = SIN(@dLat / 2) * SIN(@dLat / 2) + COS(RADIANS(@lat1))
                        * COS(RADIANS(@lat2)) * SIN(@dLon / 2) *SIN(@dLon / 2);
      SET @c = 2 * ASIN(MIN(SQRT(@a)));
      SET @d = @R * @c;
      RETURN @d;
END
GO

Pentultimate line of formula is junk. Is: SET @c = 2 * ASIN(MIN(SQRT(@a))); Should be: SET @c = 2 * atn2(sqrt(@a), sqrt(1-@a)); Not sure what R is here, use 6371 for approx Km differences.

Not sure about that. I believe this returns Statute Miles based on the Radius @R value

C#

using System;
namespace HaversineFormula
{
    /// <summary>
    /// The distance type to return the results in.
    /// </summary>
    public enum DistanceType { Miles, Kilometers };
    /// <summary>
    /// Specifies a Latitude / Longitude point.
    /// </summary>
    public struct Position
    {
        public double Latitude;
        public double Longitude;
    }
    class Haversine
    {
        /// <summary>
        /// Returns the distance in miles or kilometers of any two
        /// latitude / longitude points.
        /// </summary>
        /// <param name=”pos1″></param>
        /// <param name=”pos2″></param>
        /// <param name=”type”></param>
        /// <returns></returns>
        public double Distance(Position pos1, Position pos2, DistanceType type)
        {
            double R = (type == DistanceType.Miles) ? 3960 : 6371;
            double dLat = this.toRadian(pos2.Latitude - pos1.Latitude);
            double dLon = this.toRadian(pos2.Longitude - pos1.Longitude);
            double a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) +
                Math.Cos(this.toRadian(pos1.Latitude)) * Math.Cos(this.toRadian(pos2.Latitude)) *
                Math.Sin(dLon / 2) * Math.Sin(dLon / 2);
            double c = 2 * Math.Asin(Math.Min(1, Math.Sqrt(a)));
            double d = R * c;
            return d;
        }
        /// <summary>
        /// Convert to Radians.
        /// </summary>
        /// <param name="val"></param>
        /// <returns></returns>
        private double toRadian(double val)
        {
            return (Math.PI / 180) * val;
        }
    }
}

To use this code you will need 2 latitude/longitude points. I created a struct to hold the latitude/longitude points. Once we have the points, we create the Haversine class and call the Distance method, passing in the points and an enum specifying whether to return the results in Kilometers or Miles. We end up with the following:

Position pos1 = new Position();
pos1.Latitude = 40.7486;
pos1.Longitude = -73.9864;
Position pos2 = new Position();
pos2.Latitude = 24.7486;
pos2.Longitude = -72.9864;
Haversine hv = new Haversine();
double result = hv.Distance(pos1, pos2, DistanceType.Kilometers);

Excel

Public Function getDistance(latitude1, longitude1, latitude2, longitude2)
earth_radius = 6371
Pi = 3.14159265
deg2rad = Pi / 180

dLat = deg2rad * (latitude2 - latitude1)
dLon = deg2rad * (longitude2 - longitude1)

a = Sin(dLat / 2) * Sin(dLat / 2) + Cos(deg2rad * latitude1) * Cos(deg2rad * latitude2) * Sin(dLon / 2) * Sin(dLon / 2)
c = 2 * WorksheetFunction.Asin(Sqr(a))

d = earth_radius * c

getDistance = d

End Function

PHP

function getDistance($latitude1, $longitude1, $latitude2, $longitude2) {
	$earth_radius = 6371;
	
	$dLat = deg2rad($latitude2 - $latitude1);
	$dLon = deg2rad($longitude2 - $longitude1);
	
	$a = sin($dLat/2) * sin($dLat/2) + cos(deg2rad($latitude1)) * cos(deg2rad($latitude2)) * sin($dLon/2) * sin($dLon/2);
	$c = 2 * asin(sqrt($a));
	$d = $earth_radius * $c;
	
	return $d;
}

Smalltalk

lat1Rad := latitude * Float pi / 180.
lon1Rad := longitude * Float pi / 180.
lat2Rad := latitude2 * Float pi / 180.
lon2Rad := longitude2 * Float pi / 180.
earthRadius := 6371.00.
			
"Haversines"	
dLat := lat2Rad - lat1Rad.
dLon := lon2Rad - lon1Rad.
dLatSinSqrd := (dLat / 2) sin squared.
dLonSinSqrd := (dLon / 2) sin squared.
cosLatLat := lat2Rad cos * lat1Rad cos.
a := dLatSinSqrd + (cosLatLat * dLonSinSqrd).
c := 2 * a sqrt arcSin.
distance := earthRadius * c.

^distance

Tcl

proc deg2rad deg {expr {$deg*atan(1)*8/360}}

proc dist {lat1 lat2 lon1 lon2} {
 set R 6371
 set dLat [deg2rad [expr {($lat2 - $lat1)}]]
 set dLon [deg2rad [expr {($lon2 - $lon1)}]]
 set sdlat2 [expr {sin($dLat/2)}]
 set sdlon2 [expr {sin($dLon/2)}]
 set a [expr {$sdlat2*$sdlat2 + cos([deg2rad $lat1])*cos([deg2rad $lat2])*$sdlon2*$sdlon2}]
 set d [expr {2*$R*asin(sqrt($a))}]
}

NCL

undef("haversine_dist")
function haversine_dist( 
  lat1[*]: numeric, 
  lon1[*]: numeric, 
  lat2[*]: numeric, 
  lon2[*]: numeric 
  )
local a, c, delPhi, delPsi, R, pi, deg2rad
begin
; Earth radius in km, and other constants
   R = 6371.
   pi = 4.*atan(1.)  
   deg2rad = pi/180.
; check that all arrays are of the same length
  if ( \
    dimsizes(lat1) .ne. dimsizes(lat2) .or. \
    dimsizes(lon1) .ne. dimsizes(lon2) .or. \
    dimsizes(lat1) .ne. dimsizes(lon1) \
  ) then
    print("Function haversine_dist: Dimensions of arguments must agree. Exiting")
    exit
  end if
; Implement Haversine distance formula
; see http://www.codecodex.com/wiki/Calculate_Distance_Between_Two_Points_on_a_Globe 
  delPhi = (lat1 - lat2) * deg2rad
  delPsi = (lon1 - lon2) * deg2rad
  a = sin(delPhi/2.)*sin(delPhi/2.) + \\
    cos(lat1*deg2rad)*cos(lat2*deg2rad)*sin(delPsi/2.)*sin(delPsi/2.)
  c = 2.*atan2(sqrt(a), sqrt(1-a))
  dist = R * c 
  return(dist)
end