Calculate distance between two points on a globe

From CodeCodex

Revision as of 10:23, 5 January 2012 by 68.45.16.42 (Talk)

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;
   }
}

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);
}

Felt so hopeless looking for answers to my qusetions...until now.

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( lat ) ) ) AS DISTANCE  

DISTANCE in Miles:

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

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.

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

=2*R*ASIN(MIN(1,SQRT(SIN((RADIANS(lat2)-RADIANS(lat1))/2)^2+SIN((RADIANS(lon2-lon1))/2)^2*COS(RADIANS(lat1))*COS(RADIANS(lat2)))))

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))}]
}