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

Last change on this file since 5996 was 5942, checked in by ojw, 12 years ago

replace those horrible "degrees" (especially the negative ones) with proper compass-points when displaying POI lists

File size: 2.6 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
58def compassPoint(x):
59  points = ('N','NNE','NE','ENE','E','ESE','SE','SSE','S','SSW','SW','WSW','W','WNW','NW','NNW')
60  part = 360.0 / len(points)
61  x = x + 0.5 * part
62  while(x < 0):
63    x = x + 360
64  while(x > 360):
65    x = x - 360
66  point = int(x / part)
67  return(points[point])
68
69if(__name__ == "__main__"):
70  if(0):
71    # Test compass-point names
72    x = -180
73    while(x < 360):
74      print "%05.1f: %s" % (x, compassPoint(x))
75      x = x + 15
76 
77  a = (51.478,-0.4856)
78  b = (51.477,-0.4328)
79
80  print bearing(a,b)
81 
82  for h in (True,False):
83    start = clock()
84    it = 5000
85    for i in range(it):
86      d = distance(a,b,h)
87    t = 1000*(clock() - start) / it
88    print "%s: %1.3fms / iteration, gives %f" % (h and "haversine" or "normal", t, d)
Note: See TracBrowser for help on using the repository browser.