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

Last change on this file since 6424 was 6424, checked in by ole, 15 years ago

Experimented with scale factors in redfearn

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    # Meridian distance
135    A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
136    A2 = 3.0/8*(e2+e4/4+15*e6/128)
137    A4 = 15.0/256*(e4+3*e6/4)
138    A6 = 35*e6/3072
139   
140    term1 = a*A0*phi
141    term2 = -a*A2*sin2phi
142    term3 = a*A4*sin4phi
143    term4 = -a*A6*sin6phi
144
145    m = term1 + term2 + term3 + term4 #OK
146
147    if zone is not None and central_meridian is not None:
148        msg = 'You specified both zone and central_meridian. Provide only one of them'
149        raise Exception, msg
150   
151    # Zone
152    if zone is None:
153        zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
154
155    # Central meridian
156    if central_meridian is None:
157        central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
158    else:
159        zone = -1
160
161    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
162    omega2 = omega*omega
163    omega3 = omega*omega2
164    omega4 = omega2*omega2
165    omega5 = omega*omega4
166    omega6 = omega3*omega3
167    omega7 = omega*omega6
168    omega8 = omega4*omega4
169     
170    # Northing
171    term1 = nu*sinphi*cosphi*omega2/2 
172    term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24
173    term3 = nu*sinphi*cosphi5*\
174            (8*psi4*(11-24*t2)-28*psi3*(1-6*t2)+\
175             psi2*(1-32*t2)-psi*2*t2+t4-t2)*omega6/720
176    term4 = nu*sinphi*cosphi7*(1385-3111*t2+543*t4-t6)*omega8/40320
177    northing = false_northing + K0*(m + term1 + term2 + term3 + term4)
178
179    # Easting
180    term1 = nu*omega*cosphi
181    term2 = nu*cosphi3*(psi-t2)*omega3/6
182    term3 = nu*cosphi5*(4*psi3*(1-6*t2)+psi2*(1+8*t2)-2*psi*t2+t4)*omega5/120
183    term4 = nu*cosphi7*(61-479*t2+179*t4-t6)*omega7/5040
184    easting = false_easting + K0*(term1 + term2 + term3 + term4)
185   
186    return zone, easting, northing
187
188
189
190def convert_from_latlon_to_utm(points=None,
191                               latitudes=None,
192                               longitudes=None,
193                               false_easting=None,
194                               false_northing=None):
195    """Convert latitude and longitude data to UTM as a list of coordinates.
196
197
198    Input
199
200    points: list of points given in decimal degrees (latitude, longitude) or
201    latitudes: list of latitudes   and
202    longitudes: list of longitudes
203    false_easting (optional)
204    false_northing (optional)
205
206    Output
207
208    points: List of converted points
209    zone:   Common UTM zone for converted points
210
211
212    Notes
213
214    Assume the false_easting and false_northing are the same for each list.
215    If points end up in different UTM zones, an ANUGAerror is thrown.   
216    """
217
218    old_geo = Geo_reference()   
219    utm_points = []
220    if points == None:
221        assert len(latitudes) == len(longitudes)
222        points =  map(None, latitudes, longitudes)
223       
224    for point in points:
225       
226        zone, easting, northing = redfearn(float(point[0]),
227                                           float(point[1]),
228                                           false_easting=false_easting,
229                                           false_northing=false_northing)
230        new_geo = Geo_reference(zone)
231        old_geo.reconcile_zones(new_geo)       
232        utm_points.append([easting, northing])
233
234    return utm_points, old_geo.get_zone()
Note: See TracBrowser for help on using the repository browser.