IDL 读取葵花8(Himawari-8)HSD数据
发布日期:2021-07-25 13:05:08 浏览次数:19 分类:技术文章

本文共 4450 字,大约阅读时间需要 14 分钟。

可以实现直接输入.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,outdataEND;读取HSD数据并定标为albedo/BTFunction 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

转载地址:https://blog.csdn.net/qq_33339770/article/details/105252383 如侵犯您的版权,请留言回复原文章的地址,我们会给您删除此文章,给您带来不便请您谅解!

上一篇:Python 下载哨兵2影像
下一篇:IDL处理葵花8Himawari-8标准HSD数据——制作大气校正数据集(太阳角度数据集)

发表评论

最新留言

感谢大佬
[***.8.128.20]2024年03月18日 15时37分11秒