1 | """Implementation of Redfearn's formula to compute UTM projections from latitude and longitude |
---|
2 | |
---|
3 | """ |
---|
4 | from anuga.coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE |
---|
5 | from Numeric import array |
---|
6 | |
---|
7 | def 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 | |
---|
18 | def 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 | |
---|
34 | def 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) |
---|
166 | def 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 | |
---|
179 | def 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 |
---|