source: anuga_core/source/anuga/coordinate_transforms/lat_long_UTM_conversion.py @ 5421

Last change on this file since 5421 was 4450, checked in by duncan, 17 years ago

comments

File size: 9.7 KB
Line 
1#!/usr/bin/env python
2
3# Lat Long - UTM, UTM - Lat Long conversions
4#
5# see http://www.pygps.org
6#
7
8from math import pi, sin, cos, tan, sqrt
9
10#LatLong- UTM conversion..h
11#definitions for lat/long to UTM and UTM to lat/lng conversions
12#include <string.h>
13
14_deg2rad = pi / 180.0
15_rad2deg = 180.0 / pi
16
17_EquatorialRadius = 2
18_eccentricitySquared = 3
19
20_ellipsoid = [
21#  id, Ellipsoid name, Equatorial Radius, square of eccentricity
22# first once is a placeholder only, To allow array indices to match id numbers
23        [ -1, "Placeholder", 0, 0],
24        [ 1, "Airy", 6377563, 0.00667054],
25        [ 2, "Australian National", 6378160, 0.006694542],
26        [ 3, "Bessel 1841", 6377397, 0.006674372],
27        [ 4, "Bessel 1841 (Nambia] ", 6377484, 0.006674372],
28        [ 5, "Clarke 1866", 6378206, 0.006768658],
29        [ 6, "Clarke 1880", 6378249, 0.006803511],
30        [ 7, "Everest", 6377276, 0.006637847],
31        [ 8, "Fischer 1960 (Mercury] ", 6378166, 0.006693422],
32        [ 9, "Fischer 1968", 6378150, 0.006693422],
33        [ 10, "GRS 1967", 6378160, 0.006694605],
34        [ 11, "GRS 1980", 6378137, 0.00669438],
35        [ 12, "Helmert 1906", 6378200, 0.006693422],
36        [ 13, "Hough", 6378270, 0.00672267],
37        [ 14, "International", 6378388, 0.00672267],
38        [ 15, "Krassovsky", 6378245, 0.006693422],
39        [ 16, "Modified Airy", 6377340, 0.00667054],
40        [ 17, "Modified Everest", 6377304, 0.006637847],
41        [ 18, "Modified Fischer 1960", 6378155, 0.006693422],
42        [ 19, "South American 1969", 6378160, 0.006694542],
43        [ 20, "WGS 60", 6378165, 0.006693422],
44        [ 21, "WGS 66", 6378145, 0.006694542],
45        [ 22, "WGS-72", 6378135, 0.006694318],
46        [ 23, "WGS-84", 6378137, 0.00669438]
47]
48
49#Reference ellipsoids derived from Peter H. Dana's website-
50#http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.html
51#Department of Geography, University of Texas at Austin
52#Internet: pdana@mail.utexas.edu
53#3/22/95
54
55#Source
56#Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System
57#1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency
58
59#def LLtoUTM(int ReferenceEllipsoid, const double Lat, const double Long,
60#                        double &UTMNorthing, double &UTMEasting, char* UTMZone)
61
62def LLtoUTM( Lat, Long, ReferenceEllipsoid=23):
63    """
64    converts lat/long to UTM coords.  Equations from USGS Bulletin 1532
65    East Longitudes are positive, West longitudes are negative.
66    North latitudes are positive, South latitudes are negative
67    Lat and Long are in decimal degrees
68    Written by Chuck Gantz- chuck.gantz@globalstar.com
69    """
70    a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
71    eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
72    k0 = 0.9996
73
74#Make sure the longitude is between -180.00 .. 179.9
75    LongTemp = (Long+180)-int((Long+180)/360)*360-180 # -180.00 .. 179.9
76
77    LatRad = Lat*_deg2rad
78    LongRad = LongTemp*_deg2rad
79
80    ZoneNumber = int((LongTemp + 180)/6) + 1
81
82    if Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0:
83        ZoneNumber = 32
84
85    # Special zones for Svalbard
86    if Lat >= 72.0 and Lat < 84.0:
87        if  LongTemp >= 0.0  and LongTemp <  9.0:ZoneNumber = 31
88        elif LongTemp >= 9.0  and LongTemp < 21.0: ZoneNumber = 33
89        elif LongTemp >= 21.0 and LongTemp < 33.0: ZoneNumber = 35
90        elif LongTemp >= 33.0 and LongTemp < 42.0: ZoneNumber = 37
91
92    LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 #+3 puts origin in middle of zone
93    LongOriginRad = LongOrigin * _deg2rad
94
95    #compute the UTM Zone from the latitude and longitude
96    UTMZone = "%d%c" % (ZoneNumber, _UTMLetterDesignator(Lat))
97
98    eccPrimeSquared = (eccSquared)/(1-eccSquared)
99    N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad))
100    T = tan(LatRad)*tan(LatRad)
101    C = eccPrimeSquared*cos(LatRad)*cos(LatRad)
102    A = cos(LatRad)*(LongRad-LongOriginRad)
103
104    M = a*((1
105            - eccSquared/4
106            - 3*eccSquared*eccSquared/64
107            - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
108           - (3*eccSquared/8
109              + 3*eccSquared*eccSquared/32
110              + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
111           + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
112           - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad))
113
114    UTMEasting = (k0*N*(A+(1-T+C)*A*A*A/6
115                        + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
116                  + 500000.0)
117
118    UTMNorthing = (k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
119                                        + (61
120                                           -58*T
121                                           +T*T
122                                           +600*C
123                                           -330*eccPrimeSquared)*A*A*A*A*A*A/720)))
124
125    if Lat < 0:
126        UTMNorthing = UTMNorthing + 10000000.0; #10000000 meter offset for southern hemisphere
127    #UTMZone was originally returned here.  I don't know what the
128    #letter at the end was for.
129    #print "UTMZone", UTMZone
130    return (ZoneNumber, UTMEasting, UTMNorthing)
131
132
133def _UTMLetterDesignator(Lat):
134#This routine determines the correct UTM letter designator for the given latitude
135#returns 'Z' if latitude is outside the UTM limits of 84N to 80S
136#Written by Chuck Gantz- chuck.gantz@globalstar.com
137
138    if 84 >= Lat >= 72: return 'X'
139    elif 72 > Lat >= 64: return 'W'
140    elif 64 > Lat >= 56: return 'V'
141    elif 56 > Lat >= 48: return 'U'
142    elif 48 > Lat >= 40: return 'T'
143    elif 40 > Lat >= 32: return 'S'
144    elif 32 > Lat >= 24: return 'R'
145    elif 24 > Lat >= 16: return 'Q'
146    elif 16 > Lat >= 8: return 'P'
147    elif  8 > Lat >= 0: return 'N'
148    elif  0 > Lat >= -8: return 'M'
149    elif -8> Lat >= -16: return 'L'
150    elif -16 > Lat >= -24: return 'K'
151    elif -24 > Lat >= -32: return 'J'
152    elif -32 > Lat >= -40: return 'H'
153    elif -40 > Lat >= -48: return 'G'
154    elif -48 > Lat >= -56: return 'F'
155    elif -56 > Lat >= -64: return 'E'
156    elif -64 > Lat >= -72: return 'D'
157    elif -72 > Lat >= -80: return 'C'
158    else: return 'Z'    # if the Latitude is outside the UTM limits
159
160#void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone,
161#                         double& Lat,  double& Long )
162
163def UTMtoLL(northing, easting, zone, isSouthernHemisphere=True,
164            ReferenceEllipsoid=23):
165    """
166    converts UTM coords to lat/long.  Equations from USGS Bulletin 1532
167    East Longitudes are positive, West longitudes are negative.
168    North latitudes are positive, South latitudes are negative
169    Lat and Long are in decimal degrees.
170    Written by Chuck Gantz- chuck.gantz@globalstar.com
171    Converted to Python by Russ Nelson <nelson@crynwr.com>
172
173    FIXME: This is set up to work for the Southern Hemisphere.
174
175Using
176http://www.ga.gov.au/geodesy/datums/redfearn_geo_to_grid.jsp
177
178    Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
179Zone:   36   
180Easting:  511669.521  Northing: 19328195.112
181Latitude:   84  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
182Grid Convergence:  0  -59 ' 40.28 ''  Point Scale: 0.99960166
183
184____________
185Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
186Zone:   36   
187Easting:  519384.803  Northing: 1118247.585
188Latitude:   -80  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
189Grid Convergence:  0  59 ' 5.32 ''  Point Scale: 0.99960459
190___________
191Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
192Zone:   36   
193Easting:  611263.812  Northing: 10110547.106
194Latitude:   1  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
195Grid Convergence:  0  -1 ' 2.84 ''  Point Scale: 0.99975325
196______________
197Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
198Zone:   36   
199Easting:  611263.812  Northing: 9889452.894
200Latitude:   -1  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
201Grid Convergence:  0  1 ' 2.84 ''  Point Scale: 0.99975325
202
203So this uses a false northing of 10000000 in the both hemispheres.
204ArcGIS used a false northing of 0 in the northern hem though.
205Therefore it is difficult to actually know what hemisphere you are in.
206    """
207    k0 = 0.9996
208    a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
209    eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
210    e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared))
211
212    x = easting - 500000.0 #remove 500,000 meter offset for longitude
213    y = northing
214
215    ZoneNumber = int(zone)
216    if isSouthernHemisphere:
217        y -= 10000000.0         # remove 10,000,000 meter offset used
218                                # for southern hemisphere
219
220    LongOrigin = (ZoneNumber - 1)*6 - 180 + 3  # +3 puts origin in middle of zone
221
222    eccPrimeSquared = (eccSquared)/(1-eccSquared)
223
224    M = y / k0
225    mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256))
226
227    phi1Rad = (mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
228               + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
229               +(151*e1*e1*e1/96)*sin(6*mu))
230    phi1 = phi1Rad*_rad2deg;
231
232    N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad))
233    T1 = tan(phi1Rad)*tan(phi1Rad)
234    C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad)
235    R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5)
236    D = x/(N1*k0)
237
238    Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
239                                          +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720)
240    Lat = Lat * _rad2deg
241
242    Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
243            *D*D*D*D*D/120)/cos(phi1Rad)
244    Long = LongOrigin + Long * _rad2deg
245    return (Lat, Long)
246
247if __name__ == '__main__':
248    (z, e, n) = LLtoUTM(45.00, -75.00, 23)
249    print z, e, n
250    (lat, lon) = UTMtoLL(n, e, z, 23)
251    print lat, lon
Note: See TracBrowser for help on using the repository browser.