source: anuga_core/source/anuga/coordinate_transforms/redfearn.py @ 3614

Last change on this file since 3614 was 3614, checked in by ole, 18 years ago

Updated and tested new convert_from_latlon_to_utm

File size: 6.5 KB
Line 
1"""Implementation of Redfearn's formula to compute UTM projections from latitude and longitude
2
3"""
4from anuga.coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE
5from Numeric import array
6
7def degminsec2decimal_degrees(dd,mm,ss):
8    assert abs(mm) == mm
9    assert abs(ss) == ss   
10
11    if dd < 0:
12        sign = -1
13    else:
14        sign = 1
15   
16    return sign * (abs(dd) + mm/60. + ss/3600.)     
17   
18def decimal_degrees2degminsec(dec):
19
20    if dec < 0:
21        sign = -1
22    else:
23        sign = 1
24
25    dec = abs(dec)   
26    dd = int(dec)
27    f = dec-dd
28
29    mm = int(f*60)
30    ss = (f*60-mm)*60
31   
32    return sign*dd, mm, ss
33
34def redfearn(lat, lon, false_easting=None, false_northing=None):
35    """Compute UTM projection using Redfearn's formula
36
37    lat, lon is latitude and longitude in decimal degrees
38
39    If false easting and northing are specified they will override
40    the standard
41    """
42
43
44    from math import pi, sqrt, sin, cos, tan
45   
46
47
48    #GDA Specifications
49    a = 6378137.0                       #Semi major axis
50    inverse_flattening = 298.257222101  #1/f
51    K0 = 0.9996                         #Central scale factor   
52    zone_width = 6                      #Degrees
53
54    longitude_of_central_meridian_zone0 = -183   
55    longitude_of_western_edge_zone0 = -186
56
57    if false_easting is None:
58        false_easting = 500000
59
60    if false_northing is None:
61        if lat < 0:
62            false_northing = 10000000  #Southern hemisphere
63        else:
64            false_northing = 0         #Northern hemisphere)
65       
66   
67    #Derived constants
68    f = 1.0/inverse_flattening
69    b = a*(1-f)       #Semi minor axis
70
71    e2 = 2*f - f*f#    = f*(2-f) = (a^2-b^2/a^2   #Eccentricity
72    e = sqrt(e2)
73    e2_ = e2/(1-e2)   # = (a^2-b^2)/b^2 #Second eccentricity
74    e_ = sqrt(e2_)
75    e4 = e2*e2
76    e6 = e2*e4
77
78    #Foot point latitude
79    n = (a-b)/(a+b) #Same as e2 - why ?
80    n2 = n*n
81    n3 = n*n2
82    n4 = n2*n2
83
84    G = a*(1-n)*(1-n2)*(1+9*n2/4+225*n4/64)*pi/180
85
86
87    phi = lat*pi/180     #Convert latitude to radians
88
89    sinphi = sin(phi)   
90    sin2phi = sin(2*phi)
91    sin4phi = sin(4*phi)
92    sin6phi = sin(6*phi)
93
94    cosphi = cos(phi)
95    cosphi2 = cosphi*cosphi
96    cosphi3 = cosphi*cosphi2
97    cosphi4 = cosphi2*cosphi2
98    cosphi5 = cosphi*cosphi4   
99    cosphi6 = cosphi2*cosphi4
100    cosphi7 = cosphi*cosphi6
101    cosphi8 = cosphi4*cosphi4       
102
103    t = tan(phi)
104    t2 = t*t
105    t4 = t2*t2
106    t6 = t2*t4
107   
108    #Radius of Curvature
109    rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5
110    nu = a/(1-e2*sinphi*sinphi)**0.5
111    psi = nu/rho
112    psi2 = psi*psi
113    psi3 = psi*psi2
114    psi4 = psi2*psi2
115
116
117
118    #Meridian distance
119
120    A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
121    A2 = 3.0/8*(e2+e4/4+15*e6/128)
122    A4 = 15.0/256*(e4+3*e6/4)
123    A6 = 35*e6/3072
124   
125    term1 = a*A0*phi
126    term2 = -a*A2*sin2phi
127    term3 = a*A4*sin4phi
128    term4 = -a*A6*sin6phi
129
130    m = term1 + term2 + term3 + term4 #OK
131
132    #Zone
133    zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
134    central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
135
136    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
137    omega2 = omega*omega
138    omega3 = omega*omega2
139    omega4 = omega2*omega2
140    omega5 = omega*omega4
141    omega6 = omega3*omega3
142    omega7 = omega*omega6
143    omega8 = omega4*omega4
144   
145   
146    #Northing
147    term1 = nu*sinphi*cosphi*omega2/2 
148    term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24
149    term3 = nu*sinphi*cosphi5*\
150            (8*psi4*(11-24*t2)-28*psi3*(1-6*t2)+\
151             psi2*(1-32*t2)-psi*2*t2+t4-t2)*omega6/720
152    term4 = nu*sinphi*cosphi7*(1385-3111*t2+543*t4-t6)*omega8/40320
153    northing = false_northing + K0*(m + term1 + term2 + term3 + term4)
154
155    #Easting
156    term1 = nu*omega*cosphi
157    term2 = nu*cosphi3*(psi-t2)*omega3/6
158    term3 = nu*cosphi5*(4*psi3*(1-6*t2)+psi2*(1+8*t2)-2*psi*t2+t4)*omega5/120
159    term4 = nu*cosphi7*(61-479*t2+179*t4-t6)*omega7/5040
160    easting = false_easting + K0*(term1 + term2 + term3 + term4)
161   
162    return zone, easting, northing
163
164
165
166def convert_points_from_latlon_to_utm(points,
167                                      false_easting=None,
168                                      false_northing=None):
169    """Convert a list of points given in latitude and longitude to UTM
170
171
172    Input
173
174    points: list of points given in decimal degrees (latitude, longitude)
175    false_easting (optional)
176    false_northing (optional)
177
178    Output
179
180    zone:   UTM zone for converted points
181    points: List of converted points
182
183
184    Notes
185
186    Assume the false_easting and false_northing are the same for each list.
187    If points end up belonging to different UTM zones, an ANUGAerror is thrown.
188
189   
190   
191    """
192
193    old_geo = Geo_reference()   
194    utm_points = []
195    for point in points:
196        zone, easting, northing = redfearn(float(point[0]),
197                                           float(point[1]),
198                                           false_easting=false_easting,
199                                           false_northing=false_northing)
200
201        new_geo = Geo_reference(zone)
202        old_geo.reconcile_zones(new_geo)
203       
204        utm_points.append([easting, northing])
205
206
207    return old_geo.get_zone(), utm_points
208
209
210   
211
212
213# FIXME (Ole): Is this not supersedede by the above?
214def convert_lats_longs(latitudes,
215                       longitudes,
216                       false_easting=None,
217                       false_northing=None):
218    """
219    Redfearn a list of latitudes and longitudes.
220
221    Assume the false_easting and false_northing are the same for each list.
222
223    Note, if a zone turns out to be different, an ANUGAerror is thrown.
224    """
225    # Using geo_ref so there is only one point in the code where zones
226    # are checked.
227    old_geo = Geo_reference()
228    points = []
229    for lat, long in map(None, latitudes, longitudes,):
230        zone, easting, northing = redfearn(float(lat),
231                                           float(long),
232                                           false_easting=false_easting,
233                                           false_northing=false_northing)
234        new_geo = Geo_reference(zone)
235        old_geo.reconcile_zones(new_geo)
236        #print "old_geo.get_zone()", old_geo.get_zone()
237        points.append([easting, northing])
238       
239    return old_geo.get_zone(), points
Note: See TracBrowser for help on using the repository browser.