[2945] | 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) |
---|
[3069] | 63 | if (n_elements(NumEntries) eq 1) then ByteEntries = intarr(strnam, NumEntries) |
---|
[2945] | 64 | ByteId = intarr(strwmo, NumObs) |
---|
| 65 | ByteType = intarr(strtyp, NumObs) |
---|
| 66 | ByteJulDRef = intarr(strjuld) |
---|
| 67 | VarNames = strarr(NumVars) |
---|
[3069] | 68 | if (n_elements(NumEntries) eq 1) then Entries = strarr(NumEntries) |
---|
[2945] | 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') |
---|
[3069] | 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 |
---|
[2945] | 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 |
---|