read_cdfobs.pro 7.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221
  1. ;+---------------------------------------------------------------------------
  2. PRO read_cdfobs, Files, NumObs=NumObs, $
  3. Latitudes=Latitudes, Longitudes=Longitudes, Depths=Depths, $
  4. Obs=Observations, ModelVals=ModelVals, qcs=QCs, $
  5. Ob2=Observations2, ModelVal2=ModelVals2, qc2=QCs2, $
  6. Ob3=Observations3, $
  7. Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, kobs=kobs, $
  8. Salinity=Salinity, nodates=nodates, types=types, $
  9. filetype=ObsType, quiet=quiet, MDT=MDT, $
  10. ProfileNum=ProfileNum, error=error, $
  11. notfussy=notfussy,VarName=VarName
  12. ;+--------------------------------------------------------------------------
  13. ; Read in observation and feedback files
  14. ; detects filetype and calls the appropriate reading routine
  15. ;
  16. ; Author: D. J. Lea Feb 2008
  17. ;+--------------------------------------------------------------------------
  18. ; Declare error label.
  19. ON_IOERROR, IOERROR
  20. error=0
  21. title=''
  22. ; Set types to undefined
  23. types=-1 & tempvar=size(temporary(types))
  24. ; Read in netcdf data file and observation operator
  25. ;2. Work out which type of data is in the files by looking at the first one.
  26. ncid = ncdf_open(Files(0), /nowrite)
  27. if ncid lt 0 then message, 'Error opening file '+File
  28. result = NCDF_ATTINQ( ncid, 'title', /global)
  29. if result.datatype ne 'UNKNOWN' then ncdf_attget, ncid, 'title', Title, /global
  30. if string(Title) eq "NEMO observation operator output" then ObsType='feedback' $
  31. else $
  32. if string(Title) eq "Forecast class 4 file" then ObsType='ForecastClass4' $
  33. else $
  34. ObsType = 'none'
  35. if ObsType eq 'none' then begin
  36. varid = ncdf_varid(ncid, 'POTM_CORRECTED')
  37. if varid ge 0 then ObsType = 'Prof'
  38. varid = ncdf_varid(ncid, 'SLA')
  39. if varid ge 0 then ObsType = 'SSH'
  40. varid = ncdf_varid(ncid, 'SST')
  41. varid2 = ncdf_varid(ncid, 'sea_surface_temperature')
  42. varid=max([varid,varid2])
  43. if varid ge 0 then ObsType = 'SST'
  44. varid = ncdf_varid(ncid, 'SEAICE')
  45. varid2 = ncdf_varid(ncid, 'sea_ice_concentration')
  46. varid=max([varid,varid2])
  47. if varid ge 0 then ObsType = 'SEAICE'
  48. varid = ncdf_varid(ncid, 'LOGCHL')
  49. varid2 = ncdf_varid(ncid, 'LogChl')
  50. varid3 = ncdf_varid(ncid, 'CHL1_mean')
  51. varid=max([varid,varid2,varid3])
  52. if varid ge 0 then ObsType = 'LOGCHL'
  53. varid = ncdf_varid(ncid, 'PCO2')
  54. if varid ge 0 then ObsType = 'PCO2'
  55. varid = ncdf_varid(ncid, 'UCRT')
  56. varid2 = ncdf_varid(ncid, 'VCRT')
  57. varid=max([varid,varid2])
  58. if varid ge 0 then ObsType = 'CRT'
  59. endif
  60. ncdf_close, ncid
  61. if not keyword_set(quiet) then print, 'Reading in data as type ',ObsType
  62. if ObsType eq 'feedback' then begin
  63. if n_elements(DepRange) eq 0 then DepRange=[0,5000]
  64. read_feedback, Files, DepRange=DepRange, VarName=VarName, NumData=numobs, $
  65. OutLats=Latitudes, OutLons=Longitudes, $
  66. OutInstruments=Instruments, OutPlatform=Platform, $
  67. OutDeps=Depths, $
  68. OutObs1=Observations, OutObs2=Observations2, $
  69. OutMod1=ModelVals, OutMod2=ModelVals2, $
  70. OutQC1=QCs, OutQC2=QCs2, $
  71. MDT=MDT, OutDates=Dates, $
  72. rmdi=rmdi, quiet=quiet, $
  73. OutProfileNum=ProfileNum
  74. ; info, instruments
  75. ; info, platform
  76. types=Platform
  77. if n_elements(instruments) gt 0 then types=Platform+' '+Instruments
  78. print, 'Varname: ',VarName
  79. endif else if ObsType eq 'Prof' then begin
  80. if keyword_set(PlotModel) then begin
  81. OutTMod = 1
  82. OutSMod = 1
  83. endif
  84. if n_elements(DepRange) eq 0 then DepRange=[0,5000]
  85. read_enact, Files, DepRange=DepRange, NumData=NumData, $
  86. OutLats=Latitudes, OutLons=Longitudes, $
  87. Instruments=Instruments, Platform=Platform, $
  88. OutDeps=Depths, OutSObs=OutSObs, OutSQC=OutSQC, $
  89. OutTObs=OutTObs, OutTQC=OutTQC, OutDates=Dates, $
  90. OutTMod=OutTMod, OutSMod=OutSMod, rmdi=rmdi, quiet=quiet, $
  91. iobs=iobs, jobs=jobs, kobs=kobs, $
  92. ProfileNum=ProfileNum
  93. types=Instruments+' '+Platform
  94. ; if salinity keyword is set then put salinity values first
  95. if keyword_set(Salinity) then begin
  96. Observations = OutSObs
  97. ModelVals = OutSMod
  98. QCs = OutSQC
  99. Observations2 = OutTObs
  100. ModelVals2 = OutTMod
  101. QCs2 = OutTQC
  102. endif else begin
  103. Observations = OutTObs
  104. ModelVals = OutTMod
  105. QCs = OutTQC
  106. Observations2 = OutSObs
  107. ModelVals2 = OutSMod
  108. QCs2 = OutSQC
  109. endelse
  110. numobs=numdata
  111. numobs=n_elements(latitudes)
  112. endif else if ObsType eq 'SSH' then begin
  113. read_sla, Files, NumObs=NumObs, $
  114. Latitudes=Latitudes, Longitudes=Longitudes, $
  115. ObsSLA=Observations, ModSLA=ModelVals, SLAQC=QCs, $
  116. MDT=MDT, Satellites=Satellites, types=types, Dates=Dates, rmdi=rmdi, $
  117. iobs=iobs, jobs=jobs, mstp=mstp, quiet=quiet
  118. Depths = fltarr(NumObs)
  119. endif else if ObsType eq 'SST' then begin
  120. read_sst, Files, NumObs=NumObs, $
  121. Latitudes=Latitudes, Longitudes=Longitudes, $
  122. ObsSST=Observations, ModSST=ModelVals, SSTQC=QCs, $
  123. Dates=Dates, rmdi=rmdi, types=types, iobs=iobs, jobs=jobs, $
  124. quiet=quiet
  125. Depths = fltarr(NumObs)
  126. endif else if ObsType eq 'SEAICE' then begin
  127. read_seaice, Files, NumObs=NumObs, $
  128. Latitudes=Latitudes, Longitudes=Longitudes, $
  129. Obs=Observations, Modarr=ModelVals, QC=QCs, $
  130. Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
  131. quiet=quiet
  132. Depths = fltarr(NumObs)
  133. endif else if ObsType eq 'LOGCHL' then begin
  134. read_chl, Files, NumObs=NumObs, $
  135. Latitudes=Latitudes, Longitudes=Longitudes, $
  136. Obs=Observations, Modarr=ModelVals, QC=QCs, $
  137. Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
  138. quiet=quiet
  139. Depths = fltarr(NumObs)
  140. endif else if ObsType eq 'PCO2' then begin
  141. read_pco2, Files, NumObs=NumObs, $
  142. Latitudes=Latitudes, Longitudes=Longitudes, $
  143. Obs=Observations, Modarr=ModelVals, QC=QCs, $
  144. Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
  145. quiet=quiet
  146. Depths = fltarr(NumObs)
  147. endif else if ObsType eq 'CRT' then begin
  148. read_crt, Files, NumObs=NumObs, $
  149. OutLats=Latitudes, OutLons=Longitudes, $
  150. OutTypes=Types, OutDates=Dates, $
  151. OutUObs=Observations, OutVObs=Observations2, $
  152. OutUMod=ModelVals, OutVMod=ModelVals2, Quiet=Quiet, $
  153. FloatNum=FloatNum, QC=QCs, OutSpeed=Observations3
  154. rmdi=0.0
  155. QCs2=QCs
  156. Depths = fltarr(NumObs)
  157. endif else if ObsType eq 'ForecastClass4' then begin
  158. print, 'reading forecast class 4 files'
  159. read_forc_class4, Files, NumObs=NumObs, $
  160. Latitudes=Latitudes, Longitudes=Longitudes, $
  161. Depths=Depths, $
  162. Types=types, Dates=Dates, $
  163. Obsarr=Observations, Obs2arr=Observations2, $
  164. Modarr=ModelVals, Mod2arr=ModelVals2, $
  165. QCs1=QCs, QCs2=QCs2, rmdi=rmdi, $
  166. notfussy=notfussy
  167. endif else message, 'Error: ObsType is not set correctly'
  168. if (n_elements(types) eq 0) then begin
  169. types=replicate(rmdi,NumObs)
  170. endif
  171. goto, NOERROR
  172. IOERROR: error=1
  173. print,'read_cdfobs: an error occurred trying read files: ',files
  174. NOERROR:
  175. END