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 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 . """ 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.

Michal Niklas 14 years, 11 months ago

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) 14 years, 11 months ago

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

Gabriel Genellina 14 years, 11 months ago

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

a 14 years ago

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

James Dyson 13 years, 2 months ago

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 13 years, 2 months ago

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 8 years, 8 months ago

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)

### Required Modules

• (none specified)