source: subversion/applications/routing/pyroute-dev/geometry.py @ 18454

Last change on this file since 18454 was 18454, checked in by buerste, 10 years ago

-further updates of spaces to tabs

  • Property svn:keywords set to Rev
File size: 2.3 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4"""Handles geometry on the earth's surface (e.g. bearing/distance)
5
6Usage:
7        (library code)
8"""
9
10__version__ = "$Rev: 18454 $"
11__license__ = """This program is free software: you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation, either version 3 of the License, or
14(at your option) any later version.
15
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19GNU General Public License for more details.
20
21You should have received a copy of the GNU General Public License
22along with this program. If not, see <http://www.gnu.org/licenses/>."""
23_debug = 0
24
25
26from math import *
27from time import clock
28
29def bearing(a,b):
30        dlat = radians(b[0] - a[0])
31        dlon = radians(b[1] - a[1])
32
33        dlon = dlon * cos(radians(a[0]))
34
35        return(degrees(atan2(dlon, dlat)))
36
37def distance(a,b, haversine=False):
38        if(haversine):
39                return(distance_haversine(a,b))
40        lat1 = radians(a[0])
41        lon1 = radians(a[1])
42        lat2 = radians(b[0])
43        lon2 = radians(b[1])
44        d = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1)) * 6371;
45        return(d)
46
47
48def distance_haversine(a,b):
49        dlat = radians(b[0] - a[0])
50        dlon = radians(b[1] - a[1])
51        a1 = sin(0.5 * dlat)
52        a2a = cos(radians(a[0]))
53        a2b = cos(radians(b[0]))
54        a3 = sin(0.5 * dlon)
55
56        a = a1*a1 + a2a*a2b * a3*a3
57        c = 2 * atan2(sqrt(a), sqrt(1-a))
58        d = 6371 * c
59        return(d)
60
61def compassPoint(x):
62        points = ('N','NNE','NE','ENE','E','ESE','SE','SSE','S','SSW','SW','WSW','W','WNW','NW','NNW')
63        part = 360.0 / len(points)
64        x = x + 0.5 * part
65        while(x < 0):
66                x = x + 360
67        while(x > 360):
68                x = x - 360
69        point = int(x / part)
70        return(points[point])
71
72if(__name__ == "__main__"):
73        if(0):
74                # Test compass-point names
75                x = -180
76                while(x < 360):
77                        print "%05.1f: %s" % (x, compassPoint(x))
78                        x = x + 15
79
80        a = (51.478,-0.4856)
81        b = (51.477,-0.4328)
82
83        print bearing(a,b)
84
85        for h in (True,False):
86                start = clock()
87                it = 5000
88                for i in range(it):
89                        d = distance(a,b,h)
90                t = 1000*(clock() - start) / it
91                print "%s: %1.3fms / iteration, gives %f" % (h and "haversine" or "normal", t, d)
Note: See TracBrowser for help on using the repository browser.