source: subversion/applications/routing/pyroute/geometry.py @ 5890

Last change on this file since 5890 was 5742, checked in by ojw, 12 years ago

add benchmarking for distance functions

File size: 2.1 KB
Line 
1#!/usr/bin/python
2#-----------------------------------------------------------------------------
3# Handles geometry on the earth's surface (e.g. bearing/distance)
4#
5# Usage:
6#   (library code)
7#-----------------------------------------------------------------------------
8# Copyright 2007, Oliver White
9#
10# This program is free software: you can redistribute it and/or modify
11# it under the terms of the GNU General Public License as published by
12# the Free Software Foundation, either version 3 of the License, or
13# (at your option) any later version.
14#
15# This program is distributed in the hope that it will be useful,
16# but WITHOUT ANY WARRANTY; without even the implied warranty of
17# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18# GNU General Public License for more details.
19#
20# You should have received a copy of the GNU General Public License
21# along with this program.  If not, see <http://www.gnu.org/licenses/>.
22#-----------------------------------------------------------------------------
23from math import *
24from time import clock
25
26def bearing(a,b):
27  dlat = radians(b[0] - a[0])
28  dlon = radians(b[1] - a[1])
29
30  dlon = dlon * cos(radians(a[0]))
31 
32  return(degrees(atan2(dlon, dlat)))
33
34def distance(a,b, haversine=False):
35  if(haversine):
36    return(distance_haversine(a,b))
37  lat1 = radians(a[0])
38  lon1 = radians(a[1])
39  lat2 = radians(b[0])
40  lon2 = radians(b[1])
41  d = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1)) * 6371;
42  return(d)
43
44       
45def distance_haversine(a,b):
46  dlat = radians(b[0] - a[0])
47  dlon = radians(b[1] - a[1])
48  a1 = sin(0.5 * dlat)
49  a2a = cos(radians(a[0]))
50  a2b = cos(radians(b[0]))
51  a3 = sin(0.5 * dlon)
52 
53  a = a1*a1 + a2a*a2b * a3*a3
54  c = 2 * atan2(sqrt(a), sqrt(1-a))
55  d = 6371 * c
56  return(d)
57
58 
59if(__name__ == "__main__"):
60  a = (51.478,-0.4856)
61  b = (51.477,-0.4328)
62
63  print bearing(a,b)
64 
65  for h in (True,False):
66    start = clock()
67    it = 5000
68    for i in range(it):
69      d = distance(a,b,h)
70    t = 1000*(clock() - start) / it
71    print "%s: %1.3fms / iteration, gives %f" % (h and "haversine" or "normal", t, d)
Note: See TracBrowser for help on using the repository browser.