# source:anuga_core/source/anuga/coordinate_transforms/redfearn.py@5054

Last change on this file since 5054 was 5054, checked in by ole, 16 years ago

Comment about origin of redfearn implementation

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