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

Last change on this file since 3739 was 3739, checked in by duncan, 17 years ago

resolved convert lat longs

File size: 5.7 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_from_latlon_to_utm(latitudes=None,
167                               longitudes=None,
168                               points=None,
169                               false_easting=None,
170                               false_northing=None):
171    """Convert latitude and longitude data to UTM as a list of coordinates.
172
173
174    Input
175
176    points: list of points given in decimal degrees (latitude, longitude) or
177    latitudes: list of latitudes   and
178    longitudes: list of longitudes
179    false_easting (optional)
180    false_northing (optional)
181
182    Output
183
184    points: List of converted points
185    zone:   Common UTM zone for converted points
186
187
188    Notes
189
190    Assume the false_easting and false_northing are the same for each list.
191    If points end up belonging to different UTM zones, an ANUGAerror is thrown.   
192    """
193
194    old_geo = Geo_reference()   
195    utm_points = []
196    if points == None:
197       
198        assert len(latitudes) == len(longitudes)
199        points =  map(None, latitudes, longitudes)
200    for point in points:
201        zone, easting, northing = redfearn(float(point[0]),
202                                           float(point[1]),
203                                           false_easting=false_easting,
204                                           false_northing=false_northing)
205        new_geo = Geo_reference(zone)
206        old_geo.reconcile_zones(new_geo)       
207        utm_points.append([easting, northing])
208
209    return utm_points, old_geo.get_zone()
Note: See TracBrowser for help on using the repository browser.