本文共 5737 字，大约阅读时间需要 19 分钟。

代码解析全云盘FLDK全波段（16个）的数据并进行定标、重投影和裁剪，重投影用的GLT生成方法参考之前的博客https://blog.csdn.net/qq_33339770/article/details/102957857

PRO Himawari8 COMPILE_OPT idl2  start=systime(1)  inputfolder='***\HS_H08_20190916_0300_FLDK'  ;获取影像时间  basename=file_basename(inputfolder)  ;获取1km和2km定标后影像数据  imgdata=H8_Preprocess(inputfolder=inputfolder)   R10_data=imgdata['R10']  R20_data=imgdata['R20']  ;销毁哈希表  OBJ_DESTROY,imgdata  ;应用GTL重投影并裁剪  outR10=GLTproject(inputdata=R10_data,results=results,timename=timename)  outR20=GLTproject(inputdata=R20_data,results=results,timename=timename)  print,systime(1)-start    END;对合成的影像应用GLT几何校正并裁减成NC的范围[80,60,200,-60]Function GLTproject,inputdata=inputdata,results=results,timename=timename  COMPILE_OPT idl2  e=envi(/headless)    inputraster=ENVIRaster(inputdata,URI=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True'))  inputraster.Save  bands=inputraster.NBANDS  cols=inputraster.NCOLUMNS  if cols eq 11000 then begin    res='1000'    gltfile='***\GLT_1000M.dat'  endif else if cols eq 5500 then begin    res='2000'    gltfile='***\GLT_2000M.dat'  endif    gltraster=e.OpenRaster(gltfile)  glt_id=ENVIRastertoFID(gltraster)  ENVI_FILE_QUERY,glt_id,dims=glt_dims  ;应用GLT校正  fid=ENVIRasterToFID(inputraster)  out_name=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True')  ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT',FID=fid,GLT_DIMS=glt_dims,GLT_FID=glt_id,OUT_NAME=out_name,pos=lindgen(bands),r_fid=out_fid  ;裁剪输出  outraster=ENVIFIDtoRaster(out_fid)  SpatialRef = outraster.SPATIALREF  SpatialRef.ConvertLonLatToMap, 80.0, 60.0, MapX, MapY  SpatialRef.ConvertLonLatToMap, 200.0, -60.0, MapX2, MapY2  subraster = outraster.Subset(SPATIALREF=SpatialRef,SUB_RECT=[MapX, MapY2, MapX2, MapY])  subraster.Export,'******.dat','ENVI'       gltraster.Close  inputraster.Close  outraster.Close  subraster.Close  return,outtiffEND;--解析数据并将500m重采样至1000m,合成两个数据集--Function H8_Preprocess,inputfolder=inputfolder  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)    ;必须找到全部10个部分的数据    if N_ELEMENTS(band_files) EQ 10 THEN BEGIN      band_data=[]      seg=['01','02','03','04','05','06','07','08','09','10']      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]]        ;销毁哈希表        OBJ_DESTROY,segment      endfor      if i lt 2 then begin        R10_data[*,*,i]=band_data      endif else if i eq 2 then begin        e=envi(/headless);把B03五百米分辨率重采样        band_raster=ENVIRaster(band_data,URI=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True'))        band_raster.Save        sampleRaster=ENVIResampleRaster(band_raster, DIMENSIONS=[11000,11000]);, METHOD='Bilinear')        band_data=sampleRaster.GetData()        sampleRaster.Close        band_raster.Close        e.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      endelse    endif  endfor    outdata=hash('R10',R10_data,'R20',R20_data)  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,时间,太阳坐标  out=hash('pixels',outdata,'time',imgtime,'sun_pos',sun_pos)  return,out;data  END

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