| 1 | PRO read_feedback, Files, DepRange=DepRange, VarName=VarNames, NumData=NumData, $ |
|---|
| 2 | OutLats=OutLats, OutLons=OutLons, $ |
|---|
| 3 | OutInstruments=OutInst, OutPlatform=OutPlatform, $ |
|---|
| 4 | OutDeps=OutDeps, $ |
|---|
| 5 | OutObs1=OutObs1, OutObs2=OutObs2, $ |
|---|
| 6 | OutMod1=OutMod1, OutMod2=OutMod2, $ |
|---|
| 7 | OutQC1=OutQC1, OutQC2=OutQC2, $ |
|---|
| 8 | MDT=MDT, OutDates=OutDates, $ |
|---|
| 9 | rmdi=rmdi, quiet=quiet, $ |
|---|
| 10 | OutProfileNum=OutProfileNum |
|---|
| 11 | |
|---|
| 12 | ;------------------------------------------------------------------------------------------ |
|---|
| 13 | ;Program to read in data from a feedback format file |
|---|
| 14 | ; |
|---|
| 15 | ; Inputs: |
|---|
| 16 | ; File => File name to be read in |
|---|
| 17 | ; DepRange => (optional) Range of depths for which the data is to be extracted. |
|---|
| 18 | ; |
|---|
| 19 | ; Outputs: |
|---|
| 20 | |
|---|
| 21 | ; |
|---|
| 22 | ;Author: Matt Martin. 29/09/09 |
|---|
| 23 | ;------------------------------------------------------------------------------------------ |
|---|
| 24 | rmdi = 99999. |
|---|
| 25 | if NOT keyword_set(DepRange) then DepRange=[0.,1.e1] ;Default to top 10 metres. |
|---|
| 26 | NumFiles=n_elements(Files) |
|---|
| 27 | |
|---|
| 28 | ifile2=0 |
|---|
| 29 | numdata=0 |
|---|
| 30 | for ifile = 0, NumFiles-1 do begin |
|---|
| 31 | File = Files(ifile) |
|---|
| 32 | |
|---|
| 33 | ;1. Open the netcdf file |
|---|
| 34 | ncid = ncdf_open(File, /nowrite) |
|---|
| 35 | if ncid lt 0 and not keyword_set(quiet) then print, 'Error opening file '+File |
|---|
| 36 | if ncid ge 0 then begin |
|---|
| 37 | if not keyword_set(quiet) then print, 'Opened file '+File+' successfully' |
|---|
| 38 | |
|---|
| 39 | ;2. Get info about file |
|---|
| 40 | ncinfo = ncdf_inquire(ncid) |
|---|
| 41 | |
|---|
| 42 | ;3. Read in the relevant dimensions |
|---|
| 43 | NumExtra = 0 |
|---|
| 44 | for idim = 0, ncinfo.ndims-1 do begin |
|---|
| 45 | ncdf_diminq, ncid, idim, name, dimlen |
|---|
| 46 | if name eq 'N_OBS' then NumObs = dimlen $ |
|---|
| 47 | else if name eq 'N_LEVELS' then NumLevels = dimlen $ |
|---|
| 48 | else if name eq 'N_VARS' then NumVars = dimlen $ |
|---|
| 49 | else if name eq 'N_QCF' then NumFlags = dimlen $ |
|---|
| 50 | else if name eq 'N_ENTRIES' then NumEntries = dimlen $ |
|---|
| 51 | else if name eq 'N_EXTRA' then NumExtra = dimlen $ |
|---|
| 52 | else if name eq 'STRINGNAM' then strnam = dimlen $ |
|---|
| 53 | else if name eq 'STRINGGRID' then strgrid = dimlen $ |
|---|
| 54 | else if name eq 'STRINGWMO' then strwmo = dimlen $ |
|---|
| 55 | else if name eq 'STRINGTYP' then strtyp = dimlen $ |
|---|
| 56 | else if name eq 'STRINGJULD' then strjuld = dimlen |
|---|
| 57 | endfor |
|---|
| 58 | |
|---|
| 59 | ;4. Set up arrays |
|---|
| 60 | |
|---|
| 61 | if (NumObs gt 0) then begin |
|---|
| 62 | ByteNam = intarr(strnam, NumVars) |
|---|
| 63 | if (n_elements(NumEntries) eq 1) then ByteEntries = intarr(strnam, NumEntries) |
|---|
| 64 | ByteId = intarr(strwmo, NumObs) |
|---|
| 65 | ByteType = intarr(strtyp, NumObs) |
|---|
| 66 | ByteJulDRef = intarr(strjuld) |
|---|
| 67 | VarNames = strarr(NumVars) |
|---|
| 68 | if (n_elements(NumEntries) eq 1) then Entries = strarr(NumEntries) |
|---|
| 69 | Id = strarr(NumObs) |
|---|
| 70 | Type = strarr(NumObs) |
|---|
| 71 | if NumExtra gt 0 then begin |
|---|
| 72 | ByteExtra = intarr(strnam, NumExtra) |
|---|
| 73 | Extra = strarr(NumExtra) |
|---|
| 74 | endif |
|---|
| 75 | Longitude = fltarr(NumObs) |
|---|
| 76 | Latitude = fltarr(NumObs) |
|---|
| 77 | Depth = fltarr(NumLevels, NumObs) |
|---|
| 78 | DepQC = intarr(NumLevels, NumObs) |
|---|
| 79 | JulD = dblarr(NumObs) |
|---|
| 80 | ObsQC = intarr(NumObs) |
|---|
| 81 | PosQC = intarr(NumObs) |
|---|
| 82 | JulDQC = intarr(NumObs) |
|---|
| 83 | Obs = fltarr(NumLevels, NumObs, NumVars) |
|---|
| 84 | Hx = fltarr(NumLevels, NumObs, NumVars) |
|---|
| 85 | VarQC = intarr(NumObs, NumVars) |
|---|
| 86 | LevelQC = intarr(NumLevels, NumObs, NumVars) |
|---|
| 87 | |
|---|
| 88 | ;5. Read in the data |
|---|
| 89 | varid = ncdf_varid(ncid, 'VARIABLES') |
|---|
| 90 | ncdf_varget, ncid, varid, ByteNam |
|---|
| 91 | ; info, VarNames |
|---|
| 92 | ; info, ByteNam |
|---|
| 93 | for ivar= 0, NumVars-1 do begin |
|---|
| 94 | VarNames(ivar) = string(ByteNam(*,ivar)) |
|---|
| 95 | endfor |
|---|
| 96 | |
|---|
| 97 | if ifile2 eq 0 then VarName = VarNames(0) |
|---|
| 98 | if VarName ne VarNames(0) then message, 'Can only read in from files containing the same variables' |
|---|
| 99 | |
|---|
| 100 | varid = ncdf_varid(ncid, 'JULD_REFERENCE') |
|---|
| 101 | ncdf_varget, ncid, varid, ByteJulDRef |
|---|
| 102 | JulDRef = string(ByteJulDRef) |
|---|
| 103 | RefDate = JULDAY(fix(strmid(JulDRef,4,2)), fix(strmid(JulDRef,6,2)), fix(strmid(JulDRef,0,4)), $ |
|---|
| 104 | fix(strmid(JulDRef,8,2)), fix(strmid(JulDRef,10,2)), fix(strmid(JulDRef,12,2))) |
|---|
| 105 | print, 'RefDate: ',RefDate |
|---|
| 106 | varid = ncdf_varid(ncid, 'JULD') |
|---|
| 107 | ncdf_varget, ncid, varid, JulD |
|---|
| 108 | ; print, JulD |
|---|
| 109 | Dates = dblarr(NumLevels, NumObs) |
|---|
| 110 | for iob = 0L, long(NumObs)-1L do begin |
|---|
| 111 | dt = RefDate + JulD(iob) |
|---|
| 112 | ; print,dt, JulD(iob) |
|---|
| 113 | Dates(0:NumLevels-1, iob) = replicate(dt, NumLevels) |
|---|
| 114 | endfor |
|---|
| 115 | |
|---|
| 116 | varid = ncdf_varid(ncid, 'STATION_IDENTIFIER') |
|---|
| 117 | ncdf_varget, ncid, varid, ByteId |
|---|
| 118 | Identifier = string(ByteId) |
|---|
| 119 | |
|---|
| 120 | varid = ncdf_varid(ncid, 'STATION_TYPE') |
|---|
| 121 | ncdf_varget, ncid, varid, ByteType |
|---|
| 122 | Type = string(ByteType) |
|---|
| 123 | |
|---|
| 124 | varid = ncdf_varid(ncid, 'LONGITUDE') |
|---|
| 125 | ncdf_varget, ncid, varid, Longitude |
|---|
| 126 | |
|---|
| 127 | varid = ncdf_varid(ncid, 'LATITUDE') |
|---|
| 128 | ncdf_varget, ncid, varid, Latitude |
|---|
| 129 | |
|---|
| 130 | varid = ncdf_varid(ncid, 'DEPTH') |
|---|
| 131 | ncdf_varget, ncid, varid, Depth |
|---|
| 132 | ncdf_attget, ncid, varid, '_Fillvalue', FillVal |
|---|
| 133 | pts = where(Depth eq FillVal, count) |
|---|
| 134 | if count gt 0 then Depth(pts) = rmdi |
|---|
| 135 | |
|---|
| 136 | varid = ncdf_varid(ncid, 'OBSERVATION_QC') |
|---|
| 137 | ncdf_varget, ncid, varid, ObsQC |
|---|
| 138 | ncdf_attget, ncid, varid, '_Fillvalue', FillVal |
|---|
| 139 | pts = where(ObsQC eq FillVal, count) |
|---|
| 140 | if count gt 0 then ObsQC(pts) = rmdi |
|---|
| 141 | |
|---|
| 142 | for ivar = 0, NumVars-1 do begin |
|---|
| 143 | |
|---|
| 144 | varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_OBS') |
|---|
| 145 | ncdf_varget, ncid, varid, tmp |
|---|
| 146 | Obs(*,*,ivar) = tmp |
|---|
| 147 | ncdf_attget, ncid, varid, '_Fillvalue', FillValObs |
|---|
| 148 | |
|---|
| 149 | varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_Hx') |
|---|
| 150 | if (varid gt -1) then begin |
|---|
| 151 | ncdf_varget, ncid, varid, tmp |
|---|
| 152 | Hx(*,*,ivar) = tmp |
|---|
| 153 | ncdf_attget, ncid, varid, '_Fillvalue', FillValHx |
|---|
| 154 | endif else begin |
|---|
| 155 | FillValHx=0 |
|---|
| 156 | endelse |
|---|
| 157 | |
|---|
| 158 | varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_QC') |
|---|
| 159 | ncdf_varget, ncid, varid, tmp |
|---|
| 160 | VarQC(*,ivar) = tmp |
|---|
| 161 | ncdf_attget, ncid, varid, '_Fillvalue', FillValQC |
|---|
| 162 | |
|---|
| 163 | |
|---|
| 164 | varid = ncdf_varid(ncid, strtrim(VarNames(ivar))+'_LEVEL_QC') |
|---|
| 165 | ncdf_varget, ncid, varid, tmp |
|---|
| 166 | LevelQC(*,*,ivar) = tmp |
|---|
| 167 | ncdf_attget, ncid, varid, '_Fillvalue', FillValLevQC |
|---|
| 168 | |
|---|
| 169 | endfor |
|---|
| 170 | |
|---|
| 171 | ; print,' DJL levelqc(*,218,0) ',levelqc(*,218,0) |
|---|
| 172 | ; print,' DJL levelqc(*,218,1) ',levelqc(*,218,1) |
|---|
| 173 | |
|---|
| 174 | pts = where(Obs eq FillValObs, count) |
|---|
| 175 | if count gt 0 then Obs(pts) = rmdi |
|---|
| 176 | pts = where(Hx eq FillValHx, count) |
|---|
| 177 | if count gt 0 then Hx(pts) = rmdi |
|---|
| 178 | pts = where(VarQC eq FillValQC, count) |
|---|
| 179 | if count gt 0 then VarQC(pts) = rmdi |
|---|
| 180 | pts = where(LevelQC eq FillValLevQC, count) |
|---|
| 181 | if count gt 0 then LevelQC(pts) = rmdi |
|---|
| 182 | |
|---|
| 183 | Obs1 = Obs(*,*,0) |
|---|
| 184 | Hx1 = Hx(*,*,0) |
|---|
| 185 | VarQC1=LevelQC(*,*,0) |
|---|
| 186 | if NumVars gt 1 then begin |
|---|
| 187 | Obs2 = Obs(*,*,1) |
|---|
| 188 | Hx2 = Hx(*,*,1) |
|---|
| 189 | VarQC2=LevelQC(*,*,1) |
|---|
| 190 | endif |
|---|
| 191 | |
|---|
| 192 | if strtrim(VarName,2) eq 'SLA' then begin |
|---|
| 193 | MDT = fltarr(NumLevels, NumObs) |
|---|
| 194 | print, NumLevels, NumObs |
|---|
| 195 | varid = ncdf_varid(ncid, 'MDT') |
|---|
| 196 | ncdf_varget, ncid, varid, MDT |
|---|
| 197 | ncdf_attget, ncid, varid, '_Fillvalue', FillVal |
|---|
| 198 | pts = where(MDT eq FillVal, count) |
|---|
| 199 | if count gt 0 then MDT(pts) = rmdi |
|---|
| 200 | Obs2 = MDT |
|---|
| 201 | endif |
|---|
| 202 | |
|---|
| 203 | ;6. Put the data into the correct form to be output from this routine |
|---|
| 204 | Platformsa = strarr(NumLevels, NumObs) |
|---|
| 205 | Instrumentsa = strarr(NumLevels, NumObs) |
|---|
| 206 | Lons = fltarr(NumLevels,NumObs) |
|---|
| 207 | Lats = fltarr(NumLevels,NumObs) |
|---|
| 208 | ProfileNum = lonarr(NumLevels,NumObs) |
|---|
| 209 | |
|---|
| 210 | for i=0,NumLevels-1 do begin |
|---|
| 211 | Platformsa(i,*)=Type |
|---|
| 212 | Instrumentsa(i,*)=Identifier |
|---|
| 213 | endfor |
|---|
| 214 | |
|---|
| 215 | for i=0L,long(NumObs)-1L do begin |
|---|
| 216 | Lats(*,i)=Latitude(i) |
|---|
| 217 | Lons(*,i)=Longitude(i) |
|---|
| 218 | ProfileNum(*,i) = i+1 |
|---|
| 219 | endfor |
|---|
| 220 | ; |
|---|
| 221 | pts = where(Depth ge DepRange(0) and Depth le DepRange(1), NumDataInc) |
|---|
| 222 | NumData=NumData+NumDataInc |
|---|
| 223 | ; |
|---|
| 224 | if ifile2 eq 0 then begin |
|---|
| 225 | OutObs1 = [Obs1(pts)] |
|---|
| 226 | OutMod1 = [Hx1(pts)] |
|---|
| 227 | OutQC1 = [VarQC1(pts)] |
|---|
| 228 | OutDeps = [Depth(pts)] |
|---|
| 229 | OutLats = [Lats(pts)] |
|---|
| 230 | OutLons = [Lons(pts)] |
|---|
| 231 | OutDates= [Dates(pts)] |
|---|
| 232 | OutInst= [Instrumentsa(pts)] |
|---|
| 233 | OutPlatform= [Platformsa(pts)] |
|---|
| 234 | OutProfileNum = [ProfileNum(pts)] |
|---|
| 235 | if NumVars eq 2 then begin |
|---|
| 236 | OutObs2 = [Obs2(pts)] |
|---|
| 237 | OutMod2 = [Hx2(pts)] |
|---|
| 238 | OutQC2 = [VarQC2(pts)] |
|---|
| 239 | endif |
|---|
| 240 | endif else begin |
|---|
| 241 | OutObs1 = [OutObs1, Obs1(pts)] |
|---|
| 242 | OutMod1 = [OutMod1, Hx1(pts)] |
|---|
| 243 | OutQC1 = [OutQC1, VarQC1(pts)] |
|---|
| 244 | OutDeps = [OutDeps, Depth(pts)] |
|---|
| 245 | OutLats = [OutLats, Lats(pts)] |
|---|
| 246 | OutLons = [OutLons, Lons(pts)] |
|---|
| 247 | OutDates= [OutDates, Dates(pts)] |
|---|
| 248 | OutInst= [OutInst, Instrumentsa(pts)] |
|---|
| 249 | OutPlatform= [OutPlatform, Platformsa(pts)] |
|---|
| 250 | OutProfileNum = [OutProfileNum, ProfileNum(pts)] |
|---|
| 251 | if NumVars eq 2 then begin |
|---|
| 252 | OutObs2 = [OutObs2, Obs2(pts)] |
|---|
| 253 | OutMod2 = [OutMod2, Hx2(pts)] |
|---|
| 254 | OutQC2 = [OutQC2, VarQC2(pts)] |
|---|
| 255 | endif |
|---|
| 256 | endelse |
|---|
| 257 | |
|---|
| 258 | if NumVars eq 1 then begin |
|---|
| 259 | if VarName ne 'SLA' then OutObs2 = fltarr(n_elements(OutObs1)) + rmdi |
|---|
| 260 | OutMod2 = fltarr(n_elements(OutMod1)) + rmdi |
|---|
| 261 | OutQC2 = fltarr(n_elements(OutQC1)) + rmdi |
|---|
| 262 | endif |
|---|
| 263 | ; |
|---|
| 264 | ifile2=ifile2+1 |
|---|
| 265 | endif |
|---|
| 266 | ncdf_close, ncid |
|---|
| 267 | |
|---|
| 268 | endif |
|---|
| 269 | |
|---|
| 270 | endfor ;ifile |
|---|
| 271 | |
|---|
| 272 | pts1 = where(OutQC1 eq 1, count1) |
|---|
| 273 | pts2 = where(OutQC1 ne 1, count2) |
|---|
| 274 | if count1 gt 0 then OutQc1(pts1) = 0 |
|---|
| 275 | if count2 gt 0 then OutQc1(pts2) = 1 |
|---|
| 276 | if NumVars gt 1 then begin |
|---|
| 277 | pts1 = where(OutQC2 eq 1, count1) |
|---|
| 278 | pts2 = where(OutQC2 ne 1, count2) |
|---|
| 279 | if count1 gt 0 then OutQc2(pts1) = 0 |
|---|
| 280 | if count2 gt 0 then OutQc2(pts2) = 1 |
|---|
| 281 | endif |
|---|
| 282 | |
|---|
| 283 | END |
|---|
