-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathutm_zone.pro
72 lines (64 loc) · 2.06 KB
/
utm_zone.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
;+
; NAME:
; UTM_ZONE
;
; PURPOSE:
; The function returns the UTM zone number for a given longitude.
;
; CALLING SEQUENCE:
; result = UTM_ZONE(CM)
;
; ARGUMENTS:
; The longitude for the desired zone. Must be in decimal degrees
; and must be -180. > Lon < 180. where longitude west is negative.
; This argument may be a vector.
;
; KEYWORDS:
; CM Set this keyword equal to a named variable to return the central meridian for each
; value in the input argument.
; DOUBLE Set this keyword to force calculations in double precision.
;
; EXAMPLE:
;
; To find the UTM zone for -67 degrees (67 West)...
; IDL> Print, UTM_ZONE(-67)
; 19
;
; 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.
;
; see pages 57-58 and Figure 11
;
; MODIFICATION HISTORY:
; Written by Ben Tupper, 13JUL2000.
; 27 JAN 2003
; Added DOUBLE. BT
; Patched up problem with lon = -180 => Zone = 0
; now it forces the zone to be between 1 and 60 inclusive
; Removed loop in zone calculation.
; 2008-05-20 BT (IDL 6.3) Changed zone calulation to use VALUE_LOCATE
; and prtected the input argument (what was I thinking!).
;-
FUNCTION UTM_ZONE, iLon, CM = CM, Double = Double
On_Error, 2
Dbl = Keyword_set(Double)
Machine = Machar(Double = Dbl)
Lon = iLon
;be sure that the longitude is a bit smaller than extremes
CM = Dbl ? $
(-180.d + Machine.eps) > Lon < (180.d - Machine.eps) : $
(-180. + Machine.eps) > Lon < (180. - Machine.eps)
;make an array of bounding meridians for each zone
Bounds = Dbl ? $
Dindgen(61)*6.d - 180.d : $
Findgen(61)*6. - 180.
N = N_elements(CM)
;Zone = 1 > ((Where ( Bounds GE CM) )[0] ) < 60
Zone = 1 > (VALUE_LOCATE(bounds, cm) + 1) < 60
;if requested transform CM into its proper value
If Arg_Present(CM) Then CM = Bounds[Zone-1L]+3.
Return, Zone
END