source: branches/numpy/anuga/coordinate_transforms/point.py @ 6442

Last change on this file since 6442 was 2253, checked in by duncan, 19 years ago

creating flatter structure

File size: 3.5 KB
Line 
1"""point.py - Represents a generic point on a sphere as a Python object.
2
3   See documentation of class Point for details.
4   Ole Nielsen, ANU 2002
5"""
6
7
8from math import cos, sin, pi
9
10def acos(c):
11  """acos -  Safe inverse cosine
12 
13     Input argument c is shrunk to admissible interval
14     to avoid case where a small rounding error causes
15     a math domain error.
16  """
17  from math import acos
18
19  if c > 1: c=1
20  if c < -1: c=-1
21
22  return acos(c)
23
24
25
26
27class Point:
28    """Definition of a generic point on the sphere.
29 
30    Defines a point in terms of latitude and longitude
31    and computes distances to other points on the sphere.
32
33    Initialise as
34      Point(lat, lon), where lat and lon are in decimal degrees (dd.dddd)
35     
36
37    Public Methods:
38        distance_to(P)
39        bearing_to(P)
40        dist(P)
41
42    Author: Ole Nielsen, ANU 2002       
43    """
44     
45    # class constants
46    R = 6372000  # Approximate radius of Earth (m)
47    degrees2radians = pi/180.0
48   
49
50    def __init__(self, latitude, longitude):
51
52        assert(latitude >= -90 and latitude <= 90.0)
53        assert(longitude >= -180 and longitude <= 180.0)
54
55        self.latitude =  float(latitude)
56        self.longitude = float(longitude) 
57
58        lat = latitude * self.degrees2radians    # Converted to radians
59        lon = longitude * self.degrees2radians   # Converted to radians
60        self.coslat = cos(lat)
61        self.coslon = cos(lon)
62        self.sinlat = sin(lat)
63        self.sinlon = sin(lon)
64
65    def BearingTo(self,P):
66        """ Bearing (in degrees) to point P"""
67        AZ = self.AZ(P)
68        return int(round(AZ/self.degrees2radians))
69
70   
71    def DistanceTo(self,P):
72        """ Distance to point P"""           
73        GCA = self.GCA(P)
74        return self.R*GCA
75
76
77    def Dist(self,P):
78        """ Very cheap and rough approximation to distance"""
79        return max(abs(self.latitude-P.latitude),abs(self.longitude-P.longitude))
80
81
82    #--------------------------------------------------------------------------
83
84    def __repr__(self):
85        d = 2
86        lat = round(self.latitude,d)
87        lon = round(self.longitude,d)
88        return ' (' + str(lat)+ ', '+ str(lon)+')'
89
90
91    def GCA(self,P):
92        """ Compute the Creat Circle Angle (GCA) between current point and P.
93        """
94       
95        alpha = P.coslon*self.coslon + P.sinlon*self.sinlon
96        # The original formula is alpha = cos(self.lon - P.lon)
97        # but rewriting lets us make us of precomputed trigonometric values.
98       
99        x = alpha*self.coslat*P.coslat + self.sinlat*P.sinlat
100        return acos(x)
101
102   
103    def AZ(self,P):
104        """ Compute Azimuth bearing (AZ) from current point to P"""
105
106        # Compute cosine(AZ), where AZ is the azimuth angle
107        GCA = self.GCA(P)
108        c = P.sinlat - self.sinlat*cos(GCA)
109        c = c/self.coslat/sin(GCA)
110
111        AZ = acos(c)
112
113        # Reverse direction if bearing is westward, i.e. sin(self.lon - P.lon) > 0
114        # Without this correction the bearing due west, say, will be 90 degrees
115        # because the formulas work in the positive direction which is east.
116        #
117        # Precomputed trigonometric values are used to rewrite the formula:
118
119        if self.sinlon*P.coslon - self.coslon*P.sinlon > 0: AZ = 2*pi - AZ
120   
121        return AZ
122
123
124
125#-----------------------------------------------------------------------
126
127if __name__ == "__main__":
128    pass
Note: See TracBrowser for help on using the repository browser.