Context Navigation

read_grads.pro @ 223

Last change on this file since 223 was 223, checked in by pinsard, 17 years ago
improvements of some *.pro
Property svn:keywords set to `Id`
File size: 11.0 KB

Line
1	;------------------------------------------------------------
2	;------------------------------------------------------------
3	;------------------------------------------------------------
4	;+
5	; @file_comments
6	; reading grads file (except "data type station" or "grib")
7	; from the grads control file even if there is multiple data files.
8	;
9	; @categories
10	; Reading
11	;
12	; @param VAR {in}{required}
13	; the variable name
14	;
15	; @param DATE1 {in}{required}
16	; date of the beginning (yyyymmdd if TIMESTEP is not activate)
17	;
18	; @param DATE2 {in}{optional}
19	; last date. Optional, if not specified date2=date1
20	;
21	; @keyword FILENAME
22	; the grads control file name: 'xxxx.ctl'
23	;
24	; @keyword TIMESTEP
25	; to specify that the dates are time steps instead of true calendar.
26	;
27	;---------------
28	; NOT yet available
29	;---------------
30	;
31	; @keyword BOX {type=A 4 or 6 elements 1d array, [lon1,lon2,lat1,lat2, depth1, depth2]}
32	; It specifies the area where data must be read
33	;
34	; @keyword EVERYTHING
35	;
36	;
37	; @keyword NOSTRUCT
38	;
39	;
40	; @keyword _EXTRA
41	; Used to pass your keywords
42	;
43	; @returns
44	; an array
45	;
46	; @uses
47	; common
48	;
49	; @restrictions
50	; define all the grid parameters (defined in common.pro)
51	; associated to the data.
52	;
53	; this function call the procedure scanfile that use the
54	; unix commands grep and sed
55	;
56	; @examples
57	; IDL> a=read_grads('sst',19900101,19900131,filename='outputs.ctl')
58	; IDL> plt, a
59	;
60	; @history
61	; Sebastien Masson (smasson\@lodyc.jussieu.fr)
62	;
63	; @version
64	; $Id$
65	;
66	;-
67	;------------------------------------------------------------
68	;------------------------------------------------------------
69	;------------------------------------------------------------
70
71	FUNCTION read_grads, var, date1, date2, FILENAME = filename, BOX=box, TIMESTEP = timestep, EVERYTHING = everything, NOSTRUCT = nostruct, _EXTRA = ex
72	;---------------------------------------------------------
73	;
74	compile_opt idl2, strictarrsubs
75	;
76	@cm_4mesh
77	@cm_4data
78	@cm_4cal
79	IF NOT keyword_set(key_forgetold) THEN BEGIN
80	@updatenew
81	ENDIF
82	;------------------------------------------------------------
83	; we find the filename.
84	;------------------------------------------------------------
85	filename = isafile(FILENAME = filename, IODIRECTORY = iodir, _EXTRA = ex)
86	if size(filename, /type) NE 7 then $
87	return, report('read_ncdf cancelled')
88	;------------------------------------------------------------
89	; we scan the control file called filename
90	;------------------------------------------------------------
91	scanctl, filename, filesname, jpt1file, varsname, varslev, swapbytes, bigendian, littleendian, f77sequential, fileheader, theader, xyheader, VARFMT = varfmt, _EXTRA = ex
92	if n_elements(varfmt) EQ 0 then varfmt = 'float'
93	;------------------------
94	;------------------------
95	; check date1 and date2 and found the starting index "t1" and the
96	; ending index "t2" that corresponds to the time series specified by
97	; date1 and date2 for the time axis defined in the .ctl file.
98	;------------------------
99	;------------------------
100	if n_elements(date1) EQ 0 then begin
101	t0 = 0
102	t1 = 0
103	ENDIF
104	if n_elements(date2) EQ 0 then date2 = date1
105	if keyword_set(timestep) then BEGIN
106	if date1 GT date2 then begin
107	print, 'date2 must be larger than date1'
108	return, -1
109	endif
110	t1 = 0 > long(date1) < (jpt-1)
111	t2 = t1 > long(date2) < (jpt-1)
112	ENDIF ELSE BEGIN
113	jdate1 = date2jul(date1, /grads) < time[jpt-1]
114	jdate2 = time[0] > date2jul(date2, /grads)
115	if jdate1 GT jdate2 then begin
116	print, 'date2 must be larger than date1'
117	return, -1
118	endif
119	t1 = (where(time GE jdate1))[0]
120	tmp = where(time LE jdate2, t2)
121	t2 = t2-1
122	ENDELSE
123	if t2 LT t1 then begin
124	print, 'There is no date between date1 and date2'
125	return, -1
126	endif
127	jpt2read = t2-t1+1
128	;------------------------
129	;------------------------
130	; index of the variable
131	;------------------------
132	;------------------------
133	varid = where(strlowcase(varsname) EQ strlowcase(var))
134	varid = varid[0]
135	if varid EQ -1 then begin
136	print, var+' not found in the variable list of '+filename
137	return, -1
138	ENDIF
139	varname = var
140	if varslev[varid] EQ 1 then res = fltarr(jpi, jpj, jpt2read, /nozero) $
141	ELSE res = fltarr(jpi, jpj, varslev[varid], jpt2read, /nozero)
142	;------------------------
143	;------------------------
144	; find the first file to be read according to the file list, the
145	; number of time step in each file and t1 and t2
146	;------------------------
147	;------------------------
148	indf2read = t1/jpt1file
149	startread = t1-indf2read*jpt1file
150	alreadyread = 0
151	readagain:
152	jpt2read1file = min([jpt1file-startread, jpt2read])
153	f2read = filesname[indf2read]
154	;------------------------
155	;------------------------
156	; opening
157	;------------------------
158	;------------------------;
159	; check the existence of the file
160	f2read = isafile(filename = f2read, iodirectory = iodir, _EXTRA = ex)
161	; if the file is stored on tape
162	if !version.os_family EQ 'unix' then spawn, 'file '+f2read+' > /dev/null'
163	; open the file
164	openr, unit, f2read, /get_lun, error=err $
165	, swap_if_little_endian = bigendian $
166	, swap_if_big_endian = littleendian $
167	, swap_endian = swapbytes
168	if err ne 0 then begin
169	print,!err_string
170	return, -1
171	endif
172	;------------------------
173	case varfmt of
174	'byte':fmtsz = 1l
175	'uint':fmtsz = 2l
176	'int':fmtsz = 2l
177	'long':fmtsz = 4l
178	'float':fmtsz = 4l
179	endcase
180	;------------------------
181	; check its size
182	addf77sec = long(42f77sequential)
183	xyblocsize = xyheader + addf77sec(xyheader NE 0) + jpijpj*fmtsz +addf77sec
184	nxybloc = long(total(varslev))
185	filesize = (fileheader + addf77sec*(fileheader NE 0)) $
186	+(theader+addf77sec(theader NE 0) + nxyblocxyblocsize)*jpt1file
187	infof2read=fstat(unit)
188	if infof2read.size NE filesize then begin
189	print, 'According to '+filename+' the file size must be '+strtrim(filesize, 1)+' instead of '+strtrim(infof2read.size, 1)
190	print, 'jpi: '+strtrim(jpi, 2)
191	print, 'jpj: '+strtrim(jpj, 2)
192	print, 'jpt: '+strtrim(jpt, 2)
193	print, 'format size in byte: '+strtrim(fmtsz, 2)
194	print, 'number of xy arrays: '+strtrim(nxybloc, 2)
195	return, -1
196	endif
197	;------------------------
198	;------------------------
199	; reading
200	;------------------------
201	;------------------------;
202	;------------------------;
203	; loop on the time steps to be read in one file
204	for i = 0, jpt2read1file-1 do begin
205	; computing the offset
206	offset = (fileheader + addf77sec*(fileheader NE 0)) $
207	+(theader+addf77sec(theader NE 0) + nxyblocxyblocsize)*(startread+i) $
208	+theader+addf77sec*(theader NE 0)
209	if varid NE 0 THEN $
210	offset = offset + long(total(varslev[0:varid-1]))*xyblocsize
211	; if there is only one level
212	IF varslev[varid] EQ 1 then begin
213	case varfmt of
214	'byte':a=assoc(unit, bytarr(jpi,jpj,/nozero), offset+4*f77sequential)
215	'uint':a=assoc(unit,uintarr(jpi,jpj,/nozero), offset+4*f77sequential)
216	'int':a=assoc(unit, intarr(jpi,jpj,/nozero), offset+4*f77sequential)
217	'long':a=assoc(unit, lonarr(jpi,jpj,/nozero), offset+4*f77sequential)
218	'float':a=assoc(unit,fltarr(jpi,jpj,/nozero), offset+4*f77sequential)
219	endcase
220	res[, , i+alreadyread]=a[0]
221	ENDIF ELSE BEGIN ; more than 1 level to be read
222	if f77sequential then BEGIN ; sequential access
223	case varfmt of
224	'byte':a=assoc(unit, bytarr(jpi*jpj+8, varslev[varid],/nozero), offset)
225	'uint':a=assoc(unit,uintarr(jpi*jpj+4, varslev[varid],/nozero), offset)
226	'int':a=assoc(unit, intarr(jpi*jpj+4, varslev[varid],/nozero), offset)
227	'long':a=assoc(unit, lonarr(jpi*jpj+2, varslev[varid],/nozero), offset)
228	'float':a=assoc(unit,fltarr(jpi*jpj+2, varslev[varid],/nozero), offset)
229	endcase
230	tmp = a[0]
231	case varfmt OF ; we cut the headers and tailers of f77 write
232	'byte': tmp = tmp[4:jpijpj+3, ]
233	'uint': tmp = tmp[2:jpijpj+1, ]
234	'int': tmp = tmp[2:jpijpj+1, ]
235	'long': tmp = tmp[1:jpijpj+0, ]
236	'float':tmp = tmp[1:jpijpj+0, ]
237	endcase
238	if keyword_set(key_zreverse) then res[, , , i+alreadyread]=reverse(reform(tmp, jpi, jpj, varslev[varid], /over), 3) ELSE res[, , , i+alreadyread]=reform(tmp, jpi, jpj, varslev[varid], /over)
239	ENDIF ELSE BEGIN ; direct acces
240	case varfmt of
241	'byte':a=assoc(unit, bytarr(jpi,jpj,varslev[varid],/nozero),offset)
242	'uint':a=assoc(unit,uintarr(jpi,jpj,varslev[varid],/nozero),offset)
243	'int':a=assoc(unit, intarr(jpi,jpj,varslev[varid],/nozero),offset)
244	'long':a=assoc(unit, lonarr(jpi,jpj,varslev[varid],/nozero),offset)
245	'float':a=assoc(unit,fltarr(jpi,jpj,varslev[varid],/nozero),offset)
246	endcase
247	if keyword_set(key_zreverse) then res[, , , i+alreadyread]=reverse(a[0], 3) ELSE res[, , , i+alreadyread]=a[0]
248	ENDELSE
249	ENDELSE
250	endfor
251	;------------------------
252	; close the file
253	free_lun,unit
254	close,unit
255	;------------------------
256	;------------------------
257	; do we need to read a new file to complete the time series ?
258	;------------------------
259	;------------------------
260	if jpt2read1file NE jpt2read then BEGIN
261	indf2read = indf2read+1
262	startread = 0
263	alreadyread = alreadyread+jpt2read1file
264	jpt2read = jpt2read-jpt2read1file
265	GOTO, readagain
266	ENDIF
267	;------------------------
268	;------------------------
269	; post processing
270	;------------------------
271	;------------------------
272	if keyword_set(key_yreverse) then res = reverse(res, 2)
273	if keyword_set(key_shift) then begin
274	case (size(res))[0] of
275	2:res = shift(res, key_shift, 0)
276	3:res = shift(res, key_shift, 0, 0)
277	4:res = shift(res, key_shift, 0, 0, 0)
278	endcase
279	endif
280	;------------------------
281	; mask
282	; IF varslev[varid] EQ 1 then begin
283	; if abs(valmask) LE 1e5 then notgood = where(res[, , 0] EQ valmask) $
284	; ELSE notgood = where(abs(res[, , 0]) GE abs(valmask)/10)
285	; if notgood[0] NE -1 then tmask[notgood] = 0b
286	; ENDIF ELSE BEGIN
287	; if abs(valmask) LE 1e5 then notgood = where(res[, , *, 0] EQ valmask) $
288	; ELSE notgood = where(abs(res[, , *, 0]) GE abs(valmask)/10)
289	; if notgood[0] NE -1 then tmask[notgood] = 0b
290	; ENDELSE
291	if abs(valmask) LE 1e5 then notgood = where(res EQ valmask) $
292	ELSE notgood = where(abs(res) GE abs(valmask)/10)
293	if notgood[0] NE -1 THEN res[notgood] = !values.f_nan
294	; valmask = 1e20
295	; if abs(valmask) LE 1e5 then notgood = where(res EQ valmask) $
296	; ELSE notgood = where(abs(res) GE abs(valmask)/10)
297	; if notgood[0] NE -1 THEN res[notgood] = 1e20
298	; valmask = 1e20
299	triangles_list = triangule()
300	;------------------------
301	; subdomain extraction
302
303	;------------------------
304	; time arguments
305	;------------------------
306	time = time[t1:t2]
307	jpt = t2-t1+1
308	if keyword_set(timestep) then vardate = strtrim(time[0], 2) $
309	ELSE vardate = date2string(vairdate(time[0]))
310	;------------------------
311	@updateold
312	;------------------------
313
314	return, res
315	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format