read_seaice.pro 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178
  1. PRO read_seaice, Files, NumObs=NumObs, $
  2. Latitudes=Latitudes, Longitudes=Longitudes, $
  3. Obs=Obs, modarr=modarr, qc=qcarr, $
  4. Dates=Dates, rmdi=rmdi, iobs=outiobs, jobs=outjobs, $
  5. nodates=nodates, quiet=quiet
  6. ;------------------------------------------------------------
  7. ;IDL program to read in netcdf files of sea ice data.
  8. ;
  9. ;Author: D. J. Lea Feb 2008
  10. ;
  11. ;------------------------------------------------------------
  12. rmdi = -99999.
  13. NumFiles=n_elements(Files)
  14. ;RefDate = '1950-01-01'
  15. ;!DATE_SEPARATOR='-'
  16. RefDate=JULDAY(1,1,1950,0,0) ; should be at 0 UTC
  17. ; could read in from file
  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. ;2. Read in the dimensions in the file
  28. ncinfo = ncdf_inquire(ncid)
  29. ncdf_diminq, ncid, 0, name, NData
  30. if Ndata gt 0 then begin
  31. ;------------------------------------------------------------
  32. ;3. Allocate the data arrays and read in the data
  33. lons = dblarr(NData)
  34. lats = dblarr(NData)
  35. obsval = fltarr(NData)
  36. modval = fltarr(NData)
  37. bytQC = bytarr(NData)
  38. QC = fltarr(NData)
  39. Dats = dblarr(NData)
  40. dts = replicate(!dt_base, NData)
  41. iobs = intarr(NData)
  42. jobs = intarr(NData)
  43. ; output attribute and variable info
  44. ; ncattinfo, id
  45. varid = ncdf_varid(ncid, 'lon')
  46. ; if not keyword_set(quiet) then print,varid
  47. if (varid eq -1) then varid = ncdf_varid(ncid, 'LONGITUDE')
  48. ncdf_varget, ncid, varid, lons
  49. varid = ncdf_varid(ncid, 'lat')
  50. ; if not keyword_set(quiet) then print,varid
  51. if (varid eq -1) then varid = ncdf_varid(ncid, 'LATITUDE')
  52. ncdf_varget, ncid, varid, lats
  53. if (keyword_set(nodates) eq 0) then begin
  54. varid = ncdf_varid(ncid, 'JULD')
  55. if varid ne -1 then begin
  56. ncdf_varget, ncid, varid, Dats
  57. dts = Dats+RefDate
  58. endif else begin
  59. varid = ncdf_varid(ncid, 'time')
  60. ncdf_varget, ncid, varid, secs_from_base
  61. varid = ncdf_varid(ncid, 'SeaIce_dtime')
  62. ncdf_varget, ncid, varid, dtime
  63. ncdf_attget, ncid, varid, 'scale_factor', scale_factor
  64. if not keyword_set(quiet) then print,'dtime(0): ',dtime(0), scale_factor
  65. dtime=dtime*scale_factor
  66. ; RefDate = '1981-01-01'
  67. RefDate = JULDAY(1,1,1981,0,0) ; should be ref from 0 UTC
  68. dtime = dtime + secs_from_base
  69. dts = RefDate + dtime/86400.
  70. endelse
  71. endif
  72. if not keyword_set(quiet) then print,'reading sea ice data'
  73. varid = ncdf_varid(ncid, 'sea_ice_concentration')
  74. if not keyword_set(quiet) then print,varid
  75. if (varid eq -1) then begin
  76. varid2 = ncdf_varid(ncid, 'SEAICE')
  77. if (varid2 ne -1) then ncdf_varget, ncid, varid2, obsval
  78. endif
  79. if (varid ne -1) then ncdf_varget, ncid, varid, obsval
  80. if not keyword_set(quiet) then print,'scale_factor'
  81. scale_factor=1.
  82. if (varid ne -1) then ncdf_attget, ncid, varid, 'scale_factor', scale_factor
  83. if not keyword_set(quiet) then print,'_FillValue'
  84. FillValue=99999
  85. ; ncdf_attget, ncid, varid, '_FillValue', FillValue
  86. if not keyword_set(quiet) then print,FillValue
  87. pts = where(obsval eq FillValue, count)
  88. if not keyword_set(quiet) then print,'scale_factor ',scale_factor
  89. obsval=obsval*scale_factor
  90. if count gt 0 then obsval(pts) = rmdi
  91. if not keyword_set(quiet) then print,'reading sea ice model values'
  92. varid = ncdf_varid(ncid, 'SEAICE_Hx')
  93. if (varid ne -1) then ncdf_varget, ncid, varid, modval
  94. varid = ncdf_varid(ncid, 'IOBS')
  95. if (varid ne -1) then ncdf_varget, ncid, varid, iobs
  96. varid = ncdf_varid(ncid, 'JOBS')
  97. if (varid ne -1) then ncdf_varget, ncid, varid, jobs
  98. scale_factor=1
  99. pts = where(modval eq FillValue or modval eq -9999, count)
  100. modval=modval*scale_factor
  101. if count gt 0 then modval(pts) = rmdi
  102. if not keyword_set(quiet) then print,'sea ice qc'
  103. varid = ncdf_varid(ncid, 'SEAICE_QC')
  104. if (varid ne -1) then begin
  105. if not keyword_set(quiet) then print,'bytQC'
  106. ncdf_varget, ncid, varid, bytQC
  107. ncdf_attget, ncid, varid, '_FillValue', FillValue
  108. QC(*) = 0.
  109. pts = where(bytQC eq FillValue, count)
  110. if count gt 0 then QC(pts) = rmdi
  111. pts = where(bytQC ne 48, count)
  112. if count gt 0 then QC(pts) = 1.
  113. endif else begin
  114. varid = ncdf_varid(ncid, 'confidence_flag')
  115. ncdf_varget, ncid, varid, QC
  116. endelse
  117. if ifile2 eq 0 then begin
  118. Latitudes = [float(lats)]
  119. Longitudes = [float(lons)]
  120. Obs = [obsval]
  121. Modarr = [modval]
  122. QCarr = [QC]
  123. Dates = [dts]
  124. if (n_elements(iobs) gt 0) then outiobs = [iobs]
  125. if (n_elements(jobs) gt 0) then outjobs = [jobs]
  126. endif else begin
  127. Latitudes = [Latitudes, float(lats)]
  128. Longitudes = [Longitudes, float(lons)]
  129. Obs = [Obs, obsval]
  130. Modarr = [Modarr, modval]
  131. QCarr = [QCarr, QC]
  132. Dates = [Dates, dts]
  133. if (n_elements(iobs) gt 0) then outiobs = [outiobs, iobs]
  134. if (n_elements(jobs) gt 0) then outjobs = [outjobs, jobs]
  135. endelse
  136. ifile2=ifile2+1
  137. endif
  138. if not keyword_set(quiet) then print,'closing file'
  139. ncdf_close, ncid
  140. endif
  141. endfor ; ifile
  142. NumObs = n_elements(Latitudes)
  143. if (n_elements(Modarr) ne NumObs) then Modarr=replicate(rmdi,NumObs)
  144. END