# source:anuga_core/source/anuga/coordinate_transforms/redfearn.py@3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 5.2 KB
Line
1"""Implementation of Redfearn's formula to compute UTM projections from latitude and longitude
2
3"""
4from anuga.coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE
5from Numeric import array
6
7def 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
18def 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
34def 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
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
164def convert_lats_longs(latitudes,
165                       longitudes,
166                       false_easting=None,
167                       false_northing=None):
168    """
169    Redfearn a list of latitudes and longitudes.
170
171    Assume the false_easting and false_northing are the same for each list.
172
173    Note, if a zone turns out to be different, an ANUGAerror is thrown.
174    """
175    # Using geo_ref so there is only one point in the code where zones
176    # are checked.
177    old_geo = Geo_reference()
178    points = []
179    for lat, long in map(None, latitudes, longitudes,):
180        zone, easting, northing = redfearn(float(lat),
181                                           float(long),
182                                           false_easting=false_easting,
183                                           false_northing=false_northing)
184        new_geo = Geo_reference(zone)
185        old_geo.reconcile_zones(new_geo)
186        #print "old_geo.get_zone()", old_geo.get_zone()
187        points.append([easting, northing])
188
189    return old_geo.get_zone(), points
Note: See TracBrowser for help on using the repository browser.