-
Notifications
You must be signed in to change notification settings - Fork 3
/
utm_to_ll.pro
172 lines (155 loc) · 5.29 KB
/
utm_to_ll.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
;+
; NAME:
; UTM_TO_LL
;
; PURPOSE:
; This function converts UTM coordinates to decimal degrees of lon/lat.
; Many map data are available.
;
; REFERENCE:
; J.P. Snyder, "Map projections - A working manual',1987,
; U.S.G.S. Professional Paper 1395, Supt. of Docs No: I 19.16:1395,
; U.S. Govt Printing Office,Washington, DC 20402.
;
; CATEGORY:
; Mapping.
;
; CALLING SEQUENCE:
; Result = UTM_TO_LL(East,North, [DatumName])
;
; INPUTS:
; East/North Identically sized vectors containing the UTM coordinates
; in meters for false Easting and Northing.
;
; OPTIONAL INPUTS:
; DatumName Set this equal to a scalar string of a map datum to use.
; Possible choices are:
; NAD27
; WGS84
; NAD83
; GRS80
; WGS82
; AUS1965
; KRAS1940
; INT1924
; HAY1909
; CLARKE1880
; CLARKE1866
; AIRY1830
; BESSEL1841
; EVEREST1830
;
; If this is not set and the Datum keyword is not set then NAD27/CLARKE1866 is used.
;
; KEYWORD PARAMETERS:
; Zone Set this keyword to the Zone number for the UTM coordinate grid. If not known,
; please consult the reference above. The default is 31.
; As a quick reference, zones are 6 degrees of longitude wide starting at Zone 1
; from -180 to -174 degrees. Zone 60 extends from 174 to 180 degrees.
; NoFalse Set this keyword to indicate that Easting is not in False Easting units.
; CM Set this keyword equal to a named variable to retreive the central meridian
; in degrees for each lon/lat pair.
; P Set this keyword to a named variable to retrieve the map datum parameters.
;
; OUTPUTS:
; A 2,n elements floating point array where [0,*] = Latitude and [1,*] = Longitude in
; decimal degrees.
;
; OPTIONAL OUTPUTS:
; None.
;
; COMMON BLOCKS:
; None.
;
; SIDE EFFECTS:
; None known.
;
; RESTRICTIONS:
; Calls MAP_DATUM function.
;
; EXAMPLE:
; The following example is from the worked solutions in the Snyder reference.
; Snyder's answers are Lon = -73.5 (West) and Lat = 40.5 (North)
;
;IDL> Print, UTM_TO_LL(East,North, 'CLARKE1866', Zone = 18,CM = CM)
; -73.499995 40.499994
;
; COMMENTS:
; I suspect the difference between the solution in the procedure and that in
; Snyder's worked example lies in the calculation of Eccentricity (P.E).
; Snyder gives the formula E = SQRT(2*F - F^2)...
; Thus for The 'CLARKE1866' map datum E = 0.0822717. Surprisingly,
; the worked example uses E = 0.0822719 without explanation.
; For the CLARKE1866 map datum I could substitute the Snyder's value for
; the value resulting from the formula, but I'm not inclined to do so until
; I get confirmation from the USGS that one or the other is 'wrong'. BT 14JUL2000
;
; MODIFICATION HISTORY:
; Written by: Ben Tupper Spring 1998
; Pemaqid River Company
; email: [email protected]
; tel: (207) 563 - 1048
; 248 Lower Round Pond Road
; POB 106
; Bristol, ME 04539-0106
;
; SEP 1999 Cleaned up documentation.
; 14JULY2000 Added Zone keyword default to 31. BT
; Added MAP_DATUM fucntion to support all map data listed
; on table 1 (Ch 3) in Snyder Reference.
; Added CM keyword.
; Added NoFalse Keyword.
; Cleaned up the formatting.
; Added P keyword.
; Now returns single precision.
; Datum keyword no longer supported.
;
; 14 JAN 2008 Calculations in double precision and changed
; default DATUM to 'WGS84'. DWF
;-
FUNCTION UTM_TO_LL, East, North, DatumName,$
Zone = Zone, CM=CM, P=P, NoFalse = NoFalse
On_Error, 2
If N_Params() LT 2 Then Begin
Message, 'Two inputs arguments are required'
Return, -1
EndIf
Num = N_elements(East)
If N_elements(North) NE Num Then Begin
Message,'Input vectors must be same length'
Return, -1
EndIf
X=Double(East)
Y=Double(North)
If Not KeyWord_Set(NoFalse) Then X = X - 500000.0d
If N_elements(DatumName) EQ 0 Then $
P = Map_Datum('WGS84') Else $
P = Map_Datum(StrTrim(DatumName[0],2))
E2p = P.E^2*(1.-P.E^2)
E1 = (1. - SQRT(1. - P.E^2))/(1. + SQRT(1. - P.E^2))
K0 = 0.9996D
South = Where(North GT 10000000.d, SouthCount)
M0 = FltArr(Num)
If SouthCount GT 0 Then M0[South] = 10000000.d
M = M0 + Y/K0
u = M/(P.A*(1.0-P.E^2/4.0 - 3.0*P.E^4/64.0- 5.0*P.E^6/256.0)) ;rectifying latitude
If N_elements(Zone) EQ 0 Then Zone = Replicate(31,Num)
CM = -180.0 + 3.0d*(2.*(1>Zone<60) - 1.)
;Footprint latitude in degrees
Phi1 = u + (3.0*E1/2.0 - 27.0*E1^3/32.0)*sin(2.0*u) + $
(21.0*E1^2/16.0 - 55.0*E1^4/32.0)*sin(4.0*u) + $
(151.0*E1^3/96.0)*sin(6.0*u) + $
(1097.0*E1^4/512.0)*sin(8.0*u)
C1 = E2p*(cos(Phi1))^2
T1 = (tan(Phi1))^2
N1 = P.A/SQRT(1.0 -P.E^2*(sin(Phi1))^2 )
R1 = P.A*(1.0-P.E^2)/(1.0-P.E^2*(sin(Phi1))^2)^1.5
D = X/(N1*K0)
Y = Phi1 - $
(N1*tan(Phi1)/R1)*(D^2/2.0 - (5.0 + 3.0*T1 + 10.0*C1 -4.0*C1^2 - 9.0*E2p)*D^4/24 +$
(61.0 + 90.0*T1 + 298.0*C1 + 45.0*T1^2 + 252.0*E2p - 3.0*C1^2)*D^6/720.0)
X = CM*!DtoR + $
(D - (1.0 + 2.0*T1 + C1)*D^3/6.0 + $
(5.0 - 2.0*C1 + 28.0*T1 - 3.0*C1^2 + 8.0*E2p + 24.0*T1^2)*(D^5)/120.0)/cos(Phi1)
RETURN, Transpose([[X],[Y]])*!RaDeg
end