source: utils/tools_CPL_OASIS_ticket2379_2020/OBSTOOLS/dataplot/read_sst.pro

Last change on this file was 3026, checked in by djlea, 9 years ago

Remove unneeded sec_to_dt from dataplot. Also add obstools build instructions to documentation.

File size: 6.0 KB
Line 
1PRO  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;------------------------------------------------------------
12rmdi = -99999.
13NumFiles=n_elements(Files)
14;print,NumFiles
15;RefDate = '1950-01-01'
16;!DATE_SEPARATOR='-'
17RefDate=JULDAY(1,1,1950,0,0)    ; should be at 0 UTC
18ifile2=0
19for 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
189endfor ; ifile 
190 
191NumObs = n_elements(Latitudes)
192
193END
Note: See TracBrowser for help on using the repository browser.