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

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

Reintroduced long descriptive name

File size: 5.6 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# FIXME (Ole)
166def convert_points_from_latlon_to_utm(latitudes, longitudes,
167                                      false_easting=None,
168                                      false_northing=None):
169    """Wrapper
170    """
171   
172    return convert_lats_longs(latitudes,
173                              longitudes,
174                              false_easting,
175                              false_northing)
176
177   
178
179def convert_lats_longs(latitudes,
180                       longitudes,
181                       false_easting=None,
182                       false_northing=None):
183    """
184    Redfearn a list of latitudes and longitudes.
185
186    Assume the false_easting and false_northing are the same for each list.
187
188    Note, if a zone turns out to be different, an ANUGAerror is thrown.
189    """
190    # Using geo_ref so there is only one point in the code where zones
191    # are checked.
192    old_geo = Geo_reference()
193    points = []
194    for lat, long in map(None, latitudes, longitudes,):
195        zone, easting, northing = redfearn(float(lat),
196                                           float(long),
197                                           false_easting=false_easting,
198                                           false_northing=false_northing)
199        new_geo = Geo_reference(zone)
200        old_geo.reconcile_zones(new_geo)
201        #print "old_geo.get_zone()", old_geo.get_zone()
202        points.append([easting, northing])
203       
204    return old_geo.get_zone(), points
Note: See TracBrowser for help on using the repository browser.