Context Navigation

read_grads.pro @ 150

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

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

Download in other formats:

Original Format