Welcome, guest | Sign In | My Account | Store | Cart

Implementation of Haversine formula for calculating distance between points on a sphere. Given geographic coordinates, returns distance in kilometers.

Python, 78 lines
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#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.atan2(math.sqrt(a),  math.sqrt(1-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)

Returns an orthodromic distance, assuming the Earth is an exact sphere. Not recommended for travel planning ;) Can be useful e.g. for recalculating GPS data.

7 comments

Michal Niklas 12 years, 6 months ago  # | flag

Thanks for nice recipe. I have just search such routines to calculate distance from Google Maps .kmz files, and I used your routines in my recipe: http://code.activestate.com/recipes/576782/

One warning, according to:

http://docs.python.org/reference/lexical_analysis.html#keywords

as is a keyword and cannot be used in recalculate_coordinate() parameters list.

Bartek Górny (author) 12 years, 6 months ago  # | flag

Thanks for the warning - it worked for me because I still use Python 2.5...

Gabriel Genellina 12 years, 6 months ago  # | flag

(the usual way is to use a trailing underscore (as_) because a leading underscore (_as) means "private")

a 11 years, 7 months ago  # | flag

"c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))" can be simplified to "c = 2 * math.asin(math.sqrt(a))"

James Dyson 10 years, 10 months ago  # | flag

Sorry guys, I am such a newbie at python but not programming in general. I know nearly all the basic expressions and progressing well.

Ok from the top.... I have designed a weather balloon with a fully fledged GPS tracking system but have found a small fault in the Omni-directional antenna maths. My solution is to use a Yagi Antenna but that brings up another problem. How are you going to aim it? My second solution uses a robotic aimer to aim it.

Since I am so keen on doing it myself and also do not have the budget to buy a system like this. So with a bit of fun...do it myself.

So far so good...Everything has gone to plan so far...the maths, the design, the budget, and some of the python math programming as well.

I made a simple equation to do the vertical aim: Start

Alongitude = input("Enter longitude of Antenna GPS: ")

Alatitude = input("Enter latitude of Antenna GPS: ")

Blongitude = input("Enter longitude of Balloon GPS: ")

Blatitude = input("Enter latitude of Balloon GPS: ")

aone = Blatitude - Alatitude

bone = Blongitude - Alongitude

print aone

print bone

cone = aone * 2 + bone * 2 / 2

print cone

dddone = bone/aone

ddone = math.atan(dddone)

done = math.degrees(ddone)

Baltitude = input("Enter height of Balloon GPS: ")

Aaltitude = input("Enter height of Antenna GPS: ")

distance = cone * 2 + Baltitude * 2 / 2

Caltitude = Baltitude - Aaltitude

a = Baltitude/Caltitude

b = math.atan(a)

c = math.degrees(b)

print "The degree of vertical angle is:"

print c

print "The distance between the Balloon GPS and the Antenna GPS is:"

print distance

print "The degrees of horizontal angle is:"print ddone End

Note: the end result and the data entered will be automated

Now I just had to find or make a equation to sort out the latitude and longitude to find the distance to be used in the first equation and the bearing to aim the robot antenna.

But I can't somehow just get my head around this equation here. Funnily enough have gotten confused over the data input, the type of input (I should know by now) and the actual equation myself.

I have just turned 15 (dec) so excuse any mistakes I make or my lack of understanding but I have researched but cannot understand completely.

Could anyone explain this code to me in more detail and how to implement this to my design?

Sorry about the code it when wild and clamped up

Thanks James Dyson 15

James Dyson 10 years, 10 months ago  # | flag

Made a fatal error in my last comment, the code, this is the right code.

Start

import math

Aaltitude = input("Enter height of A: ")

Opposite = input("Enter height of B: ")

Base = input("Enter base distance between A and B: ")

distance = Base * 2 + Opposite * 2 / 2

Caltitude = Opposite - Aaltitude

a = Opposite/Base

b = math.atan(a)

c = math.degrees(b)

print "The degree of vertical angle is:"

print c

print "The distance between the Balloon GPS and the Antenna GPS is:"

print distance

End

Sorry I am such a noob...

romainbeucher 6 years, 3 months ago  # | flag

This post is rather old but as I came across an issue testing it I thought it would be good to add a precision.

The code above is valid in Python 2. but will return wrong value in Python 3 That comes from the fact that it uses the controversial "/" division operator which in python 2 returns the floor division when used with 2 integers.That is not the case in python 3 (return an approximation of the true division)...

Solution is to explicity use floor division or the "//" operator both valid in python2 and 3...

Created by Bartek Górny on Wed, 27 May 2009 (MIT)
Python recipes (4591)
Bartek Górny's recipes (1)

Required Modules

  • (none specified)

Other Information and Tasks