IDL 读取葵花8（Himawari-8）HSD数据

可以实现直接输入.DAT文件，读取获得定标完成的影像、时间、太阳位置的三个数组
inputfile为.DAT文件，例如HS_H08_20170623_0250_B01_FLDK_R10_S0110.DAT
resolution为HS_H08_20170623_0250_B01_FLDK_R后的数字，代表了分辨率
代码会自动对数据进行对应的定标（可见光波段定标为表观反射率，红外定标为亮温）
时间提前转为double型，防止后面麻烦
会将500m的波段采样到一千米，然后将数据按分辨率划分，合并

;输入文件夹（存储全部DAT文件）并进行预处理
Function H8_Preprocess,inputfolder=inputfolder,e=e
  COMPILE_OPT idl2

  R10_data=make_array(11000,11000,4,/FLOAT)
  R20_data=make_array(5500,5500,12,/FLOAT)
  
  bandnames=['01','02','03','04','05','06','07','08','09','10','11','12','13','14','15','16']
  for i=0,15 do begin
    band_files=file_search(inputfolder,'*B'+bandnames[i]+'_FLDK*.DAT',/test_regular)
    if N_ELEMENTS(band_files) EQ 10 THEN BEGIN
      band_data=[]
      obstime=[]
      sunpos=[]
      for j=0,9 do begin
        infile=band_files[j]
        segment=read_hsd(inputfile=infile[0])
        imgdata=segment['pixels']
        imgtime=segment['time']
        sun_pos=segment['sun_pos']
        band_data=[[band_data],[imgdata]]
        obstime=[[obstime],[imgtime]]
        sunpos=[[sunpos],[sun_pos]]
        ;销毁哈希表
        OBJ_DESTROY,segment
      endfor
      
      if i lt 2 then begin
        R10_data[*,*,i]=band_data
        R10_obstime=obstime
        R10_sunpos=sunpos
      endif else if i eq 2 then begin       
        band_raster=ENVIRaster(band_data,URI=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True'))
        band_raster.Save
        sampleRaster=ENVIResampleRaster(band_raster, DIMENSIONS=[11000,11000])
        band_data=sampleRaster.GetData()
        sampleRaster.Close
        band_raster.Close
        R10_data[*,*,i]=band_data
      endif else if i lt 4 then begin
        R10_data[*,*,i]=band_data
      endif else begin
        R20_data[*,*,i-4]=band_data
        R20_obstime=obstime
        R20_sunpos=sunpos
      endelse

    endif
  endfor
  outdata=hash('R10',R10_data,'R10_time',R10_obstime,'R10_sunpos',R10_sunpos,'R20',R20_data,'R20_time',R20_obstime,'R20_sunpos',R20_sunpos)
  return,outdata

END
;读取HSD数据并定标为albedo/BT
Function read_hsd,inputfile=inputfile
  COMPILE_OPT idl2
  
  spli=strsplit(file_basename(inputfile),'_')
  ;获取分辨率
  resolution=fix(strmid(file_basename(inputfile),spli[6]+1,2))
  if resolution eq 5 then begin
    cols=22000
    rows=2200
  endif else if resolution eq 10 then begin
    cols=11000
    rows=1100
  endif else if resolution eq 20 then begin
    cols=5500
    rows=550
  endif
  
  openr,h8_lun,inputfile,/get_lun
  ;获取Observation start time
  point_lun,h8_lun,46
  imgtime=dblarr(1);R8
  readu,h8_lun,imgtime
  ;获取Total header length
  point_lun,h8_lun,70
  header_length=uintarr(1);I4
  readu,h8_lun,header_length
  ;获取影像
  point_lun,h8_lun,header_length
  imgdata=uintarr(cols,rows)
  readu,h8_lun,imgdata 
  ;获取Sun's position
  point_lun,h8_lun,510
  sun_pos=dblarr(3);R8
  readu,h8_lun,sun_pos
  ;获取Band number
  point_lun,h8_lun,601
  band=uintarr(1);I2
  readu,h8_lun,band  
  ;获取Gain
  point_lun,h8_lun,617
  Gain=dblarr(1);R8
  readu,h8_lun,Gain
  ;获取Offset
  point_lun,h8_lun,625
  Offset=dblarr(1);R8
  readu,h8_lun,Offset
  ;计算radiance
  radiance=imgdata*Gain[0]+Offset[0]
  
  if band le 6 then begin;计算反射率
    ;获取radiance to albedo
    point_lun,h8_lun,633
    cc=dblarr(1);R8
    readu,h8_lun,cc
    
    ;计算反射率
    albedo=radiance*cc[0]   
    ;返回值outdata
    outdata=albedo
  endif else begin;计算亮温
    ;获取Central wave length
    point_lun,h8_lun,603
    wv=dblarr(1);R8
    readu,h8_lun,wv
    ;获取radiance to brightness temperature(c0)
    point_lun,h8_lun,633
    c0=dblarr(1);R8
    readu,h8_lun,c0
    ;获取radiance to brightness temperature(c1)
    point_lun,h8_lun,641
    c1=dblarr(1);R8
    readu,h8_lun,c1
    ;获取radiance to brightness temperature(c2)
    point_lun,h8_lun,649
    c2=dblarr(1);R8
    readu,h8_lun,c2
    ;获取Speed of light (c)
    point_lun,h8_lun,681
    c=dblarr(1);R8
    readu,h8_lun,c    
    ;获取Planck constant (h)
    point_lun,h8_lun,689
    h=dblarr(1);R8
    readu,h8_lun,h
    ;获取Boltzmann constant(k)
    point_lun,h8_lun,697
    k=dblarr(1);R8
    readu,h8_lun,k
    
    ;计算亮温
    wv=wv[0]*1e-6
    rad=radiance*1e6
    Te=h[0]*c[0]/k[0]/wv[0]/(ALOG(2*h[0]*c[0]^2/(wv[0]^5*rad)+1))
    BT=c0[0]+c1[0]*Te+c2[0]*Te^2
    
    ;返回值 outdata
    outdata=BT
  endelse
  free_lun,h8_lun 
  ;返回：albedo/BT,时间,太阳坐标
  obstime=make_array(cols,rows,value=imgtime[0]*1.0d)
  sunpos=make_array(cols,rows,3,/float)
  sunpos[*,*,0]=make_array(cols,rows,value=sun_pos[0])
  sunpos[*,*,1]=make_array(cols,rows,value=sun_pos[1])
  sunpos[*,*,2]=make_array(cols,rows,value=sun_pos[2])
  out=hash('pixels',outdata,'time',obstime,'sun_pos',sunpos)
  return,out
  
END