source: branches/numpy/anuga/coordinate_transforms/redfearn.py @ 6902

Last change on this file since 6902 was 6902, checked in by rwilson, 15 years ago

Back-merge from Numeric trunk to numpy branch.

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