IDL实现FY3D MERSI2 250M/1KM数据预处理

前言

    两个月前第一次正式使用IDL，花了几天时间才实现了MERSI2的预处理和两个温度反演算法，而且代码写的乱七八糟。然后今天参考了IDL实现风云三号D星1km数据定标及几何校正这位博主的博客，代码很整齐并且里面两个读取HDF5的函数写的很漂亮，收获很多。于是重新整理了一下之前的代码，重写了250M和1KM数据预处理部分，温度反演部分还需要修正算法之后再重写。

说明

• 结合我自己的需求只处理了Band3.4.15.18.24.25，后续获取光谱响应函数后会加上6s大气校正部分。
• 由于构建250M的GLT文件时间比较长，所以没有添加构建GLT部分，构建GLT的方法可以参考ENVI-IDL技术殿堂官方博客。
• 代码末删除临时文件的方法比较暴力，带有[envitempfile]的文件都会删除，结合实际情况，谨慎使用！

代码

运行前添加以下4个文件的路径及文件名。

  fnDat = '_0250M_MS.HDF'
  fnDatkm='_1000M_MS.HDF'
  fnLoc = '_GEOQK_MS.HDF'
  fnSz  = '_GEO1K_MS.HDF'  

pro preMERSI2_250m
  compile_opt idl2, hidden
  tic
  e = ENVI()

  fnDat = '_0250M_MS.HDF'
  fnDatkm='_1000M_MS.HDF'
  fnLoc = '_GEOQK_MS.HDF'
  fnSz  = '_GEO1K_MS.HDF'  
  
  dirname = FILE_DIRNAME(fnDat)
  date = (FILE_BASENAME(fnDat)).Extract $
    ('20[1-4][0-9][0,1][0-9][0-3][0-9]_[0-9][0-9][0-9][0-9]')
  refb3_4fn = dirname + PATH_SEP() + 'MERSI2_' + date + '_B3_4TOA.dat'
  refb15_18fn = dirname + PATH_SEP() + 'MERSI2_' + date + '_B15_18TOA.dat'
  Tbb_b24_25fn = dirname + PATH_SEP() + 'MERSI2_' + date + '_B24_25Tbb.dat'
  
  fdid = H5F_OPEN(fnDat)
    ;rawRefb1 = h5GetData(fdid,'Data/EV_250_RefSB_b1')
    ;rawRefb2 = h5GetData(fdid,'Data/EV_250_RefSB_b2')
    rawRefb3 = h5GetData(fdid,'Data/EV_250_RefSB_b3')
    rawRefb4 = h5GetData(fdid,'Data/EV_250_RefSB_b4')
    
    rawEmib24 = h5GetData(fdid,'Data/EV_250_Emissive_b24')
    rawEmib25 = h5GetData(fdid,'Data/EV_250_Emissive_b25')
    calCoef = h5GetData(fdid, '/Calibration/VIS_Cal_Coeff')
    
     dst =(h5GetAttr(fdid, 'EarthSun Distance Ratio'))[0]
    ESUN = h5GetAttr(fdid, 'Solar_Irradiance')
    tbbA = h5GetAttr(fdid, 'TBB_Trans_Coefficient_A')
    tbbB = h5GetAttr(fdid, 'TBB_Trans_Coefficient_B')
  H5F_CLOSE, fdid
  
  fdkmid = H5F_OPEN(fnDatkm)
    rawRef = h5GetData(fdkmid, '/Data/EV_1KM_RefSB')
    rawRefb15 = rawRef[*,*,10]
    rawRefb18 = rawRef[*,*,13]
  H5F_CLOSE, fdkmid
   
  ;flid = H5F_OPEN(fnLoc)
    ;lat = h5GetData(flid, 'Latitude')
    ;lon = h5GetData(flid, 'Longitude')
  ;H5F_CLOSE, flid
  
  fsid = H5F_OPEN(fnSz)
    ;SenZ = FLOAT(h5GetData(fsid, 'Geolocation/SensorZenith')*  !pi / 180 )
    SolZ = FLOAT(h5GetData(fsid, 'Geolocation/SolarZenith') * !pi / 180 )
  H5F_CLOSE, fsid
  
  SolZ_TempFn = e.GetTemporaryFilename()
  Solz_raster = e.CreateRaster(SolZ_TempFn,SolZ,INHERITS_FROM = raster)
  Solz_raster.Save
  SolZ = !null
  
  rawRefb15_TempFn = e.GetTemporaryFilename()
  rawRefb15_raster = e.CreateRaster(rawRefb15_TempFn,rawRefb15,INHERITS_FROM = raster)
  rawRefb15_raster.Save
  rawRefb15 = !null
  
  rawRefb18_TempFn = e.GetTemporaryFilename()
  rawRefb18_raster = e.CreateRaster(rawRefb18_TempFn,rawRefb18,INHERITS_FROM = raster)
  rawRefb18_raster.Save
  rawRefb18 = !null
  
  rw = N_elements(rawRefb3[0,*])
  cl = N_elements(rawRefb3[*,0])

  SolZ_r = ENVIResampleRaster(Solz_raster, DIMENSIONS=[cl,rw], METHOD='Bilinear') 
  SolZ_r = SolZ_r.GetData(bands=0)
  
  rawRefb15_r = ENVIResampleRaster(rawRefb15_raster, DIMENSIONS=[cl,rw], METHOD='Bilinear')
  rawRefb15_r = rawRefb15_r.GetData(bands=0)
  
  rawRefb18_r = ENVIResampleRaster(rawRefb18_raster, DIMENSIONS=[cl,rw], METHOD='Bilinear')
  rawRefb18_r = rawRefb18_r.GetData(bands=0)
    
  refb3_4 = MAKE_ARRAY([cl,rw, 2], type = 5)
  refb3_4[*,*,0] = (rawRefb3 * calCoef[1,2] + calCoef[0,2]) * (dst^2) / cos(SolZ_r)
  refb3_4[*,*,1] = (rawRefb4 * calCoef[1,3] + calCoef[0,3]) * (dst^2) / cos(SolZ_r)
  rawRefb3 = !null
  rawRefb4 = !null
  
  refb15_18 = MAKE_ARRAY([cl,rw, 2], type = 5)
  refb15_18[*,*,0] = (rawRefb15_r * calCoef[1,14] + calCoef[0,14]) * (dst^2) / cos(SolZ_r)
  refb15_18[*,*,1] = (rawRefb18_r * calCoef[1,17] + calCoef[0,17]) * (dst^2) / cos(SolZ_r)
  rawRefb15_r = !null
  rawRefb18_r = !null
  SolZ_r = !null
  
  wn = [933.364, 836.941]
  k = 1.191043E-5
  
  Tbb_b24_25 = MAKE_ARRAY([cl,rw, 2], type = 5)
  Tbb_b24_25[*,*,0] = tbbA[4]* $
  ((1.43878 * wn[0]) / ALOG(1 + (k*wn[0]^3) / rawEmib24)) + tbbB[4]
  Tbb_b24_25[*,*,1] = tbbA[4]* $
  ((1.43878 * wn[1]) / ALOG(1 + (k*wn[1]^3) / rawEmib25)) + tbbB[5]
  rawEmib24 = !null
  rawEmib25 = !null
  
  Solz_raster.Close
  rawRefb15_raster.Close
  rawRefb18_raster.Close
  
  refb3_4_Tempfn =  e.GetTemporaryFilename()
  refb3_4_raster = e.CreateRaster(refb3_4_Tempfn,refb3_4,INHERITS_FROM = raster)
  refb3_4_raster.Save
  refb3_4id = ENVIRasterToFID(refb3_4_raster)
  refb3_4 = !null
  
  refb15_18_Tempfn =  e.GetTemporaryFilename()
  refb15_18_raster = e.CreateRaster(refb15_18_Tempfn,refb15_18,INHERITS_FROM = raster)
  refb15_18_raster.Save
  refb15_18id = ENVIRasterToFID(refb15_18_raster)
  refb15_18 = !null
  
  Tbb_b24_25_Tempfn =  e.GetTemporaryFilename()
  Tbb_b24_25_raster = e.CreateRaster(Tbb_b24_25_Tempfn,Tbb_b24_25,INHERITS_FROM = raster)
  Tbb_b24_25_raster.Save
  Tbb_b24_25id = ENVIRasterToFID(Tbb_b24_25_raster)
  Tbb_b24_25 = !null
  
  GLTfn = ENVI_PICKFILE(title = 'select GLT file')
  ENVI_OPEN_FILE,GLTfn,R_FID = GLT_id
 
  ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT', $
  BACKGROUND = 0, $
  FID = refb3_4id, $
  GLT_FID = GLT_id, $
  OUT_NAME = refb3_4fn, $
  pos = [0:1]
  
  ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT',$
  BACKGROUND = 0, $
  FID = refb15_18id, $
  GLT_FID = GLT_id, $
  OUT_NAME = refb15_18fn, $
  pos = [0:1]
  
  ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT',$
  BACKGROUND = 0, $
  FID = Tbb_b24_25id, $
  GLT_FID = GLT_id, $
  OUT_NAME = Tbb_b24_25fn, $
  pos=[0:1]

  refb3_4_raster.Close
  refb15_18_raster.Close
  Tbb_b24_25_raster.Close
  
  refb3_4Raster = e.OpenRaster(refb3_4fn)
  refb3_4Raster.Metadata.UpdateItem, 'band names', $
    ['Band3', 'Band4']
  refb3_4Raster.Metadata.UpdateItem, 'wavelength units', 'Micrometers'
  refb3_4Raster.Metadata.AddItem, 'wavelength', $
    [0.650, 0.865]
  refb3_4Raster.WriteMetadata
  refb3_4Raster.Close
  
  refb15_18Raster = e.OpenRaster(refb15_18fn)
  refb15_18Raster.Metadata.UpdateItem, 'band names', $
    ['Band15', 'Band18']
  refb15_18Raster.Metadata.UpdateItem, 'wavelength units', 'Micrometers'
  refb15_18Raster.Metadata.AddItem, 'wavelength', $
    [0.865, 0.940]
  refb15_18Raster.WriteMetadata
  refb15_18Raster.Close
  
  Tbb_b24_25Raster = e.OpenRaster(Tbb_b24_25fn)
  Tbb_b24_25Raster.Metadata.UpdateItem, 'band names', $
    ['Band24', 'Band25']
  Tbb_b24_25Raster.Metadata.UpdateItem, 'wavelength units', 'Micrometers'
  Tbb_b24_25Raster.Metadata.AddItem, 'wavelength', $
    [10.8, 12.0]
  Tbb_b24_25Raster.WriteMetadata
  Tbb_b24_25Raster.Close
  
  CD, FILE_DIRNAME(e.GetTemporaryFilename())
  TempFiles = FILE_SEARCH('[envitempfile]*')
  n = N_elements(TempFiles)
  for i = 0,n-1 do begin
    FILE_DELETE, TempFiles[i]
  endfor
  
  toc
end

function h5GetData, fid, str
  compile_opt idl2, hidden

  str_id = H5D_OPEN(fid, str)
  slope = FLOAT(h5GetAttr(str_id, 'Slope'))
  intercept = FLOAT(h5GetAttr(str_id, 'Intercept'))
  data = FLOAT(H5D_READ(str_id))
  foreach _slope, slope, index do $
    data[*, *, index] = _slope * data[*, *, index] + intercept[index]
  H5D_CLOSE, str_id

  RETURN, data
end

function h5GetAttr, fid, str
  compile_opt idl2, hidden

  str_id = H5A_OPEN_NAME(fid, str)
  attr = H5A_READ(str_id)
  H5A_CLOSE, str_id

  RETURN, attr
end