[2945] | 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 | ;------------------------------------------------------------ |
---|
[3002] | 7 | ;IDL program to read in netcdf files of altimeter data. |
---|
[2945] | 8 | ; |
---|
[3002] | 9 | ;Author: D. J. Lea - Feb 2008 |
---|
| 10 | ; |
---|
[2945] | 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 |
---|