1 | PRO read_sst, Files, NumObs=NumObs, $ |
---|
2 | Latitudes=Latitudes, Longitudes=Longitudes, $ |
---|
3 | ObsSST=ObsSST, ModSST=ModSST, SSTQC=SSTQC, $ |
---|
4 | Dates=Dates, rmdi=rmdi, Types=Types, iobs=iobs, jobs=jobs, $ |
---|
5 | quiet=quiet |
---|
6 | ;------------------------------------------------------------ |
---|
7 | ;IDL program to read in netcdf files of altimeter data. |
---|
8 | ; |
---|
9 | ;Author: D. J. Lea - Feb 2008 |
---|
10 | ; |
---|
11 | ;------------------------------------------------------------ |
---|
12 | rmdi = -99999. |
---|
13 | NumFiles=n_elements(Files) |
---|
14 | ;print,NumFiles |
---|
15 | ;RefDate = '1950-01-01' |
---|
16 | ;!DATE_SEPARATOR='-' |
---|
17 | RefDate=JULDAY(1,1,1950,0,0) ; should be at 0 UTC |
---|
18 | ifile2=0 |
---|
19 | for ifile = 0, NumFiles-1 do begin |
---|
20 | ;------------------------------------------------------------ |
---|
21 | ;1. Open the file containing the data |
---|
22 | ncid = ncdf_open(Files(ifile), /nowrite) |
---|
23 | if ncid lt 0 and not keyword_set(quiet) then print, 'Error opening file '+Files(ifile) |
---|
24 | if ncid ge 0 then begin |
---|
25 | if not keyword_set(quiet) then print, 'Opened file '+Files(ifile)+' successfully' |
---|
26 | |
---|
27 | ;------------------------------------------------------------ |
---|
28 | ;2. Read in the dimensions in the file |
---|
29 | ncinfo = ncdf_inquire(ncid) |
---|
30 | |
---|
31 | if ncinfo.Ndims eq 1 then begin |
---|
32 | ncdf_diminq, ncid, 0, name, NData |
---|
33 | endif else if ncinfo.Ndims eq 2 then begin |
---|
34 | ncdf_diminq, ncid, 1, name, NData |
---|
35 | if (name eq "string8") then ncdf_diminq, ncid, 0, name, NData |
---|
36 | endif else if ncinfo.Ndims eq 3 then begin |
---|
37 | ncdf_diminq, ncid, 1, name, NLats |
---|
38 | ncdf_diminq, ncid, 2, name, NLons |
---|
39 | NData = NLats * NLons |
---|
40 | endif |
---|
41 | |
---|
42 | if not keyword_set(quiet) then print, 'Number of SST data points ',NData |
---|
43 | |
---|
44 | if (NData gt 0) then begin |
---|
45 | ;------------------------------------------------------------ |
---|
46 | ;3. Allocate the data arrays and read in the data |
---|
47 | lons = dblarr(NData) |
---|
48 | lats = dblarr(NData) |
---|
49 | obsval = fltarr(NData) |
---|
50 | modval = fltarr(NData) |
---|
51 | bytQC = bytarr(NData) |
---|
52 | QC = fltarr(NData) |
---|
53 | obstyp = intarr(NData) |
---|
54 | Dats = dblarr(NData) |
---|
55 | dts = replicate(!dt_base, NData) |
---|
56 | iobsa = intarr(NData) |
---|
57 | jobsa = intarr(NData) |
---|
58 | |
---|
59 | varid = ncdf_varid(ncid, 'LONGITUDE') |
---|
60 | if varid ne -1 then ncdf_varget, ncid, varid, lons $ |
---|
61 | else begin |
---|
62 | varid = ncdf_varid(ncid, 'lon') |
---|
63 | ncdf_varget, ncid, varid, lons |
---|
64 | endelse |
---|
65 | |
---|
66 | varid = ncdf_varid(ncid, 'LATITUDE') |
---|
67 | if varid ne -1 then ncdf_varget, ncid, varid, lats $ |
---|
68 | else begin |
---|
69 | varid = ncdf_varid(ncid, 'lat') |
---|
70 | ncdf_varget, ncid, varid, lats |
---|
71 | endelse |
---|
72 | |
---|
73 | varid = ncdf_varid(ncid, 'JULD') |
---|
74 | if varid ne -1 then begin |
---|
75 | ncdf_varget, ncid, varid, Dats |
---|
76 | dts=RefDate + Dats |
---|
77 | |
---|
78 | endif else begin |
---|
79 | varid = ncdf_varid(ncid, 'time') |
---|
80 | ncdf_varget, ncid, varid, secs_from_base |
---|
81 | varid = ncdf_varid(ncid, 'sst_dtime') |
---|
82 | ncdf_varget, ncid, varid, dtime |
---|
83 | ncdf_attget, ncid, varid, 'scale_factor', scale_factor |
---|
84 | dtime=dtime*scale_factor |
---|
85 | ; RefDate = '1981-01-01' |
---|
86 | RefDate = JULDAY(1,1,1981,0,0) ; should be ref from 0 UTC |
---|
87 | dtime = dtime + secs_from_base |
---|
88 | dts = RefDate + dtime/86400. |
---|
89 | endelse |
---|
90 | |
---|
91 | varid = ncdf_varid(ncid, 'SST') |
---|
92 | if varid ne -1 then begin |
---|
93 | ncdf_varget, ncid, varid, obsval |
---|
94 | obsval=float(obsval) ;ensure obsval is a floating point array |
---|
95 | ncdf_attget, ncid, varid, '_FillValue', FillValue |
---|
96 | pts = where(obsval eq FillValue, count) |
---|
97 | if count gt 0 then obsval(pts) = rmdi |
---|
98 | endif else begin |
---|
99 | varid = ncdf_varid(ncid, 'sea_surface_temperature') |
---|
100 | ncdf_varget, ncid, varid, obsval |
---|
101 | obsval=float(obsval) ;ensure obsval is a floating point array |
---|
102 | ncdf_attget, ncid, varid, '_FillValue', FillValue |
---|
103 | ncdf_attget, ncid, varid, 'add_offset', Offset |
---|
104 | ncdf_attget, ncid, varid, 'scale_factor', Scale |
---|
105 | pts = where(obsval ne FillValue, count) |
---|
106 | if count gt 0 then obsval(pts) = Offset + (obsval(pts) *Scale) |
---|
107 | pts = where(obsval eq FillValue, count) |
---|
108 | if count gt 0 then obsval(pts) = rmdi |
---|
109 | endelse |
---|
110 | |
---|
111 | varid = ncdf_varid(ncid, 'SST_Hx') |
---|
112 | if varid ne -1 then begin |
---|
113 | ncdf_varget, ncid, varid, modval |
---|
114 | ncdf_attget, ncid, varid, '_FillValue', FillValue |
---|
115 | pts = where(modval eq FillValue, count) |
---|
116 | if count gt 0 then modval(pts) = rmdi |
---|
117 | endif |
---|
118 | |
---|
119 | varid = ncdf_varid(ncid, 'SST_QC') |
---|
120 | if varid ne -1 then begin |
---|
121 | ncdf_varget, ncid, varid, bytQC |
---|
122 | ncdf_attget, ncid, varid, '_FillValue', FillValue |
---|
123 | QC(*) = 0. |
---|
124 | pts = where(bytQC eq FillValue, count) |
---|
125 | if count gt 0 then QC(pts) = rmdi |
---|
126 | pts = where(bytQC ne 48, count) |
---|
127 | if count gt 0 then QC(pts) = bytQC(pts)-48 |
---|
128 | endif else begin |
---|
129 | varid = ncdf_varid(ncid, 'confidence_flag') |
---|
130 | ncdf_varget, ncid, varid, bytQC |
---|
131 | QC = float(bytQC) |
---|
132 | endelse |
---|
133 | |
---|
134 | varid = ncdf_varid(ncid, 'SST_DATA_SOURCE') |
---|
135 | if varid ne -1 then ncdf_varget, ncid, varid, obstyp $ |
---|
136 | else begin |
---|
137 | varid = ncdf_varid(ncid, 'data_source') |
---|
138 | if varid ne -1 then ncdf_varget, ncid, varid, obstyp |
---|
139 | endelse |
---|
140 | |
---|
141 | varid = ncdf_varid(ncid, 'callsign') |
---|
142 | if varid ne -1 then begin |
---|
143 | ncdf_varget, ncid, varid, callsign |
---|
144 | tmp=strtrim(string(obstyp),2)+' '+strtrim(string(callsign),2) |
---|
145 | obstyp=tmp |
---|
146 | endif else begin |
---|
147 | varid = ncdf_varid(ncid, 'SST_CALL_SIGN') |
---|
148 | if varid ne -1 then begin |
---|
149 | ncdf_varget, ncid, varid, callsign |
---|
150 | tmp=strtrim(string(obstyp),2)+' '+strtrim(string(callsign),2) |
---|
151 | obstyp=tmp |
---|
152 | endif |
---|
153 | endelse |
---|
154 | |
---|
155 | varid = ncdf_varid(ncid, 'IOBS') |
---|
156 | if (varid ne -1) then ncdf_varget, ncid, varid, iobsa |
---|
157 | |
---|
158 | varid = ncdf_varid(ncid, 'JOBS') |
---|
159 | if (varid ne -1) then ncdf_varget, ncid, varid, jobsa |
---|
160 | |
---|
161 | if ifile2 eq 0 then begin |
---|
162 | Latitudes = [float(lats)] |
---|
163 | Longitudes = [float(lons)] |
---|
164 | ObsSST = [obsval] |
---|
165 | ModSST = [modval] |
---|
166 | SSTQC = [QC] |
---|
167 | Dates = [dts] |
---|
168 | Types = [obstyp] |
---|
169 | iobs = [iobsa] |
---|
170 | jobs = [jobsa] |
---|
171 | endif else begin |
---|
172 | Latitudes = [Latitudes, float(lats)] |
---|
173 | Longitudes = [Longitudes, float(lons)] |
---|
174 | ObsSST = [ObsSST, obsval] |
---|
175 | ModSST = [ModSST, modval] |
---|
176 | SSTQC = [SSTQC, QC] |
---|
177 | Dates = [Dates, dts] |
---|
178 | Types = [Types, obstyp] |
---|
179 | iobs = [iobs, iobsa] |
---|
180 | jobs = [jobs, jobsa] |
---|
181 | endelse |
---|
182 | |
---|
183 | ifile2=ifile2+1 |
---|
184 | endif |
---|
185 | ncdf_close, ncid |
---|
186 | |
---|
187 | endif |
---|
188 | |
---|
189 | endfor ; ifile |
---|
190 | |
---|
191 | NumObs = n_elements(Latitudes) |
---|
192 | |
---|
193 | END |
---|