Implementation of Haversine formula for calculating distance between points on a sphere. Given geographic coordinates, returns distance in kilometers.
Python, 78 lines
Download
Copy to clipboard1 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.
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.
Thanks for the warning - it worked for me because I still use Python 2.5...
(the usual way is to use a trailing underscore (as_) because a leading underscore (_as) means "private")
"c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))" can be simplified to "c = 2 * math.asin(math.sqrt(a))"
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
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...
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...