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

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

Change Numeric imports to general form - ready to change to NumPy?.

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