source: trunk/procs/nc_read.pro @ 122

Last change on this file since 122 was 121, checked in by ericg, 16 years ago

mask value now correct when scale_factor is present in netCDF file
File size: 4.9 KB
Line 
1;--------------------------- 
2; Reading any netcdf file
3;--------------------------- 
4FUNCTION nc_read, file_name, var_name, ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, NO_MEAN = no_mean
5
6; arguments = file_name, varname
7; ncdf_db=<location>:<path> or just <path>
8; mot clef = iodir time_1 time_2
9
10@common
11@com_eg
12 
13; inits
14   IF debug_w THEN print, '   '
15   IF debug_w THEN print, '   ENTER nc_read...'
16
17   IF NOT keyword_set(TIME_1) THEN time_1 =  1
18   IF NOT keyword_set(TIME_2) THEN time_2 =  time_1
19   IF debug_w THEN print, '  Before read_ncdf,  key_yreverse=', key_yreverse 
20;
21
22; decide which subdomain of data to read
23   IF debug_w THEN print, '     keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA)
24   IF keyword_set(ALL_DATA) THEN BEGIN
25      tout = 1
26   ENDIF ELSE tout = 0
27   CASE vert_type OF
28      'z' : zindex = 0
29      'level' : zindex = 1
30      ELSE: zindex = 0
31   ENDCASE
32
33; local directory
34   IF strpos(ncdf_db, ':') GE 1 THEN directory = (strsplit(ncdf_db, ':', /EXTRACT))[1] $
35    ELSE directory = ncdf_db
36
37; Find the variable's attribute
38   ncdf_getatt, directory+file_name, var_name, add_offset = add_offset, scale_factor = scale_factor, $
39    missing_value = missing_value, units = units, calendar = calendar, long_name = long_name
40
41; By default units for the Z axis are meters
42   IF n_elements(gdept) GT 1 THEN BEGIN
43      IF name_level NE '-' THEN $
44       ncdf_getatt, directory+file_name, name_level, units = units_depth ELSE units_depth = 'm'
45   ENDIF ELSE units_depth = ''
46
47; Check consistency between time_2 computed from param in cmdline and the max number
48; of time steps found in the file.
49   jpt_max = find_jptmax(file_name, IODIRECTORY = directory)
50   IF time_2-time_1+1 GT jpt_max AND jpt_max NE -1 THEN time_2 = time_1+jpt_max-1
51
52; Read the data with a call to read_ncdf
53   IF debug_w THEN print, '     '
54   IF debug_w THEN print, '     Reading raw data from ',  file_name
55
56   lec_data = read_ncdf(var_name, time_1-1, time_2-1, BOXZOOM = boxzoom, FILENAME = directory+file_name, $
57                        /TIMESTEP, TOUT = tout, /NOSTRUCT, /CONT_NOFILL, GRID = vargrid, $
58                        ZINVAR = zinvar, ZINDEX = zindex)
59   field_dim = (size(lec_data))[0]
60   IF field_dim LE 0 THEN stop,  '  Something wrong happened in read_ncdf '
61
62   IF debug_w THEN print, '     '
63   IF debug_w AND keyword_set(boxzoom) THEN print,  '    boxzoom=', boxzoom
64   IF debug_w THEN print,  '    firstxt,firstyt,firstzt=', firstxt, firstyt, firstzt
65   IF debug_w THEN print,  '    lastxt,lastyt,lastzt=', firstxt, firstyt, firstzt   
66   IF debug_w THEN print,  '    zinvar=', zinvar
67   IF debug_w THEN print, '     jpt=', jpt
68   IF debug_w THEN print, '     key_shift=', key_shift
69   IF debug_w THEN print, '     After read_ncdf, key_yreverse=',key_yreverse
70   IF debug_w THEN print, '     '
71   
72; Average data along vertical if needed and update some features
73; needed for plt (data_direc, name_suff)
74   name_suff = ''
75   data_direc = ''
76   update_data, TAB = lec_data, VNAME = var_name, BOXZOOM = boxzoom, ZUNITS = units_depth, $
77    ZINVAR = zinvar, SUFFIX = name_suff, D_DIREC = data_direc, $
78    TIME_1 = time_1, TIME_2 =  time_2, NO_MEAN = no_mean
79   field_dim = (size(lec_data))[0]
80
81; Field attributes
82   field = {name: '', data: lec_data, legend: '', units: '', origin: '', dim: 0, direc: data_direc}
83   field.name = var_name
84   field.origin = directory+file_name
85   field.legend = long_name+name_suff
86   field.units = units
87   field.dim = field_dim
88
89; replace NaN with valmask
90;   idx_t = where (~finite(field.data))
91;   IF idx_t(0) NE -1 THEN field.data(idx_t) = valmask
92
93; get valmask (might need valmask = float(string(valmask))
94   valmask = 1.e20
95   IF size(missing_value, /TYPE) EQ 4 OR size(missing_value, /TYPE) EQ 5 THEN BEGIN
96      valmask = missing_value
97      ; ensure valmask is positive
98      IF valmask LT 0 THEN BEGIN
99         print, '      *** Warning valmask is negative - changing sign: ', valmask
100         idx_t = where (field.data EQ valmask)
101         IF idx_t[0] NE -1 THEN field.data(idx_t) = -valmask
102         valmask = -valmask
103         idx_t=0                ; free memory
104      ENDIF
105   ENDIF
106
107   IF debug_w THEN print, '     valmask=',valmask
108
109; if pseudo 3d mask, read mask file and set masked points to valmask
110; still to do
111
112; min/max
113
114   chardim = ' - dim = '
115   FOR i = 1, (size(field.data))[0] DO BEGIN
116      chardim = chardim+string((size(field.data))[i], format = '(I5)')
117   ENDFOR 
118
119   idx_valmsk = (where (field.data EQ valmask))
120   idx_novalmsk = (where (field.data NE valmask))
121   minf = idx_valmsk[0] EQ -1 ? min(field.data) : min(field.data(idx_novalmsk))
122   maxf = idx_valmsk[0] EQ -1 ? max(field.data) : max(field.data(idx_novalmsk))
123
124   print, '      = ',field.legend, '    [min/max = ',minf , maxf,'  ', field.units,' - masked values = ',valmask, chardim, ']'
125
126   IF debug_w THEN print,  '  ...EXIT nc_read'
127   IF debug_w THEN print,  '  '
128   
129   return, field
130
131END
Note: See TracBrowser for help on using the repository browser.