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

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

Implemented ability to specify central meridian for non-UTM projections.
Test partially passes, but has been disabled for now.

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