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
173
174
175
176
177
178
|
# Plot global climatological temperature distribution
# from the NCEP reanalysis data.
# The data is downloaded if not found and the users wants.
######################################
### local functions ###
def nearest_index(na,val)
# returns the element of na nearest to val
# na is assumed to be 1d and monotonic
if(na[0] > na[-1]) then
reversed = true
na = na[-1..0]
else
reversed = false
end
w = (na.lt(val)).where
idx = [ (w.length==0 ? 0 : w.max), na.length-2 ].min
if ( na[idx+1]-val < val-na[idx] ) then
idx = idx+1
end
if reversed then
na = na[-1..0]
idx = na.length-1-idx
end
idx
end
#####################################
### main part ###
require "numru/dcl"
require "numru/netcdf"
include NumRu
filename = "air.mon.ltm.nc"
# < download if not found and the users wants >
if !(Dir.glob(filename)[0]) then
# file not found ==> download by ftp if the user wants
host = "ftp.cdc.noaa.gov"
path = "/Datasets/ncep.reanalysis.derived/pressure/"+filename
print "\n*** question ***\n",
" File "+filename+" is not in the current directory.\n",
" Would you like to download (ftp) it from "+host+"?\n",
" (y, n)> "
ans = gets
if ans =~ /^y/ then
print " What is your email address? (needed for anonymous ftp) > "
email = gets.chop!
require "net/ftp"
print " connecting...\n"
ftp = Net::FTP.open(host, "anonymous", email, nil)
size = ftp.size(path)
print " Size of the file is #{(size/1000)} kb. Would you really like to download?\n (y, n)> "
ans = gets
if ans =~ /^y/ then
print " now downloading...\n"
ftp.getbinaryfile(path, filename)
else
print " exit\n"
exit
end
else
print " exit\n"
exit
end
end
# < open the file and read axes >
file = NetCDF.open(filename)
var = file.var("air") # temperature
lon = file.var("lon")
lat = file.var("lat")
level = file.var("level")
time = file.var("time") # in hours
t = (time.get/720).round + 1 # --> in months
axes = [ lon, lat, level, time ]
axvals = [ lon.get, lat.get, level.get, t ]
axunits = [ lon.att("units").get.gsub("_",""),
lat.att("units").get.gsub("_",""),
level.att("units").get.gsub("_",""),
"months" ]
# < graphics >
#DCL.sglset('lbuff',false)
DCL.swlset('lwait',false)
DCL.gropn(1)
first = true
while true do
begin
## / select a 2D slice /
if (first) then
print <<-EOS
** select a slice **
List desired grid numbers of lonigutede, latitude, level, time
(set only two of them)
Example:
, , 850, 1
-- horizontal slice at 850hPa and January
135, , , 2
-- vertical slice at 135E and Feburary
EOS
DCL.grfrm
else
DCL.grfrm
print "Input next slice (C-d to quit)\n"
end
print "> "
slice = gets.chop!.split(",")
if slice.length!=4 then
raise("Slice must be 4 comma-split numerics")
end
slice.collect!{|e| # "collect!" replaces elements
if e =~ /^ *$/ then
nil
else
e.to_f
end
}
iax = []
start=[] ; last=[]
slice.each_index{|i|
if slice[i] == nil then
iax.push(i)
start.push(0) ; last.push(-1) # from the beginning to the end
else
idx = nearest_index( axvals[i], slice[i] )
start.push( idx ) ; last.push( idx )
end
}
if iax.length != 2 then
raise("Specify a 2D slice")
else
x = axvals[iax[0]]
iax[0]==2 ? xr=[x.max, x.min] : xr=[x.min, x.max]
xttl = axes[iax[0]].att("long_name").get
xunits = axunits[iax[0]]
y = axvals[iax[1]]
iax[1]==2 ? yr=[y.max, y.min] : yr=[y.min, y.max]
yttl = axes[iax[1]].att("long_name").get
yunits = axunits[iax[1]]
end
## / read the slice and plot /
v = var.get("start"=>start, "end"=>last)
shp=v.shape; shp.delete(1); v.reshape!(*shp) # delete dims of length==1
#Fig.inclpoint(x, y)
DCL.grswnd( xr[0], xr[1], yr[0], yr[1] )
DCL.grsvpt(0.2,0.9,0.2,0.9)
DCL.grstrf
DCL.ussttl(xttl," ",yttl," ")
DCL.usdaxs
DCL.uwsgxa(x)
DCL.uwsgya(y)
DCL.uelset("ltone",true)
DCL.uetone(v)
DCL.udcntz(v)
first = false
rescue
print "*Error* ", $!,"\n" # show the error message in ($!)
end
end
DCL.grcls
|
