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

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

Replaced 'print' statements with log.critical() calls.

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    return (ZoneNumber, UTMEasting, UTMNorthing)
130
131
132def _UTMLetterDesignator(Lat):
133#This routine determines the correct UTM letter designator for the given latitude
134#returns 'Z' if latitude is outside the UTM limits of 84N to 80S
135#Written by Chuck Gantz- chuck.gantz@globalstar.com
136
137    if 84 >= Lat >= 72: return 'X'
138    elif 72 > Lat >= 64: return 'W'
139    elif 64 > Lat >= 56: return 'V'
140    elif 56 > Lat >= 48: return 'U'
141    elif 48 > Lat >= 40: return 'T'
142    elif 40 > Lat >= 32: return 'S'
143    elif 32 > Lat >= 24: return 'R'
144    elif 24 > Lat >= 16: return 'Q'
145    elif 16 > Lat >= 8: return 'P'
146    elif  8 > Lat >= 0: return 'N'
147    elif  0 > Lat >= -8: return 'M'
148    elif -8> Lat >= -16: return 'L'
149    elif -16 > Lat >= -24: return 'K'
150    elif -24 > Lat >= -32: return 'J'
151    elif -32 > Lat >= -40: return 'H'
152    elif -40 > Lat >= -48: return 'G'
153    elif -48 > Lat >= -56: return 'F'
154    elif -56 > Lat >= -64: return 'E'
155    elif -64 > Lat >= -72: return 'D'
156    elif -72 > Lat >= -80: return 'C'
157    else: return 'Z'    # if the Latitude is outside the UTM limits
158
159#void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone,
160#                         double& Lat,  double& Long )
161
162def UTMtoLL(northing, easting, zone, isSouthernHemisphere=True,
163            ReferenceEllipsoid=23):
164    """
165    converts UTM coords to lat/long.  Equations from USGS Bulletin 1532
166    East Longitudes are positive, West longitudes are negative.
167    North latitudes are positive, South latitudes are negative
168    Lat and Long are in decimal degrees.
169    Written by Chuck Gantz- chuck.gantz@globalstar.com
170    Converted to Python by Russ Nelson <nelson@crynwr.com>
171
172    FIXME: This is set up to work for the Southern Hemisphere.
173
174Using
175http://www.ga.gov.au/geodesy/datums/redfearn_geo_to_grid.jsp
176
177    Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
178Zone:   36   
179Easting:  511669.521  Northing: 19328195.112
180Latitude:   84  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
181Grid Convergence:  0  -59 ' 40.28 ''  Point Scale: 0.99960166
182
183____________
184Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
185Zone:   36   
186Easting:  519384.803  Northing: 1118247.585
187Latitude:   -80  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
188Grid Convergence:  0  59 ' 5.32 ''  Point Scale: 0.99960459
189___________
190Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
191Zone:   36   
192Easting:  611263.812  Northing: 10110547.106
193Latitude:   1  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
194Grid Convergence:  0  -1 ' 2.84 ''  Point Scale: 0.99975325
195______________
196Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
197Zone:   36   
198Easting:  611263.812  Northing: 9889452.894
199Latitude:   -1  0 ' 0.00000 ''  Longitude: 34  0 ' 0.00000 ''
200Grid Convergence:  0  1 ' 2.84 ''  Point Scale: 0.99975325
201
202So this uses a false northing of 10000000 in the both hemispheres.
203ArcGIS used a false northing of 0 in the northern hem though.
204Therefore it is difficult to actually know what hemisphere you are in.
205    """
206    k0 = 0.9996
207    a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
208    eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
209    e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared))
210
211    x = easting - 500000.0 #remove 500,000 meter offset for longitude
212    y = northing
213
214    ZoneNumber = int(zone)
215    if isSouthernHemisphere:
216        y -= 10000000.0         # remove 10,000,000 meter offset used
217                                # for southern hemisphere
218
219    LongOrigin = (ZoneNumber - 1)*6 - 180 + 3  # +3 puts origin in middle of zone
220
221    eccPrimeSquared = (eccSquared)/(1-eccSquared)
222
223    M = y / k0
224    mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256))
225
226    phi1Rad = (mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
227               + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
228               +(151*e1*e1*e1/96)*sin(6*mu))
229    phi1 = phi1Rad*_rad2deg;
230
231    N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad))
232    T1 = tan(phi1Rad)*tan(phi1Rad)
233    C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad)
234    R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5)
235    D = x/(N1*k0)
236
237    Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
238                                          +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720)
239    Lat = Lat * _rad2deg
240
241    Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
242            *D*D*D*D*D/120)/cos(phi1Rad)
243    Long = LongOrigin + Long * _rad2deg
244    return (Lat, Long)
245
246if __name__ == '__main__':
247    (z, e, n) = LLtoUTM(45.00, -75.00, 23)
248    print z, e, n
249    (lat, lon) = UTMtoLL(n, e, z, 23)
250    print lat, lon
Note: See TracBrowser for help on using the repository browser.