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
