批量波段提取融合

;代码可实现多波段遥感影像到有效波段的融合提取，可消除较多波段造成的数据冗余

PRO LAYER_STACKING
  COMPILE_OPT idl2
  e = ENVI(/headless)
  ENVI, /restore_base_save_files
  ENVI_BATCH_INIT

  inpath = DIALOG_PICKFILE(title='选择输入目录',/DIRECTORY);
  outpath=DIALOG_PICKFILE(title='选择输出目录',/DIRECTORY)
  CD,inpath
  filenames = FILE_SEARCH('*.dat')
  n = N_ELEMENTS(filenames)
  PRINT, n

  PRINT, '开始处理数据 : ', SYSTIME()
  FOR i=0,n-1 DO BEGIN
    file=inpath+filenames[i]
    ENVI_OPEN_FILE, File ,r_fid=fid
    ENVI_FILE_QUERY,fid,dims=dims,ns=ns,nL=nL,nb=nb,data_type=data_type,interleave=interleave,wl=wl
    map_info=ENVI_GET_MAP_INFO(fid=fid)
    o_pro=ENVI_GET_PROJECTION(fid=fid)

    ;对检索出来的数据做波段提取融合
    layerstacking=UINTARR(ns,nL,4);根据数据类型创建输出数据类型，并设定提取波段的数量为4
    ;由于博主处理的数据时landsat8数据，hdr头文件打开时data type对应的值12指的是16位无符号整数，故这里创建时采用的UINTARR()数组创建函数
    grid_file=inpath+filenames[i]
    PRINT,grid_file
    OPENR,lun,grid_file,/GET_LUN
    arr=UINTARR(ns,nL,nb);与上述创建数组的类型保持一致
    help,arr
    READU,lun,arr
    FREE_LUN,lun
    ;需要融合哪些波段照着以下格式书写：左边的0，1,2,3指的输出文件的第1,2,3,4个波段，右边见下面注释
    layerstacking[*,*,0]=arr[*,*,1];融合源文件中的第2个波段
    layerstacking[*,*,1]=arr[*,*,2];融合源文件中的第3个波段
    layerstacking[*,*,2]=arr[*,*,3];融合源文件中的第4个波段
    layerstacking[*,*,3]=arr[*,*,4];融合源文件中的第5个波段
    wl=[wl[1],wl[2],wl[3],wl[4]];设定波长单位，1234对应上述源文件第几个波段，融了几个波段，大[]里就写几个，用,隔开
    
    outname=outpath+filenames[i]
    ENVI_WRITE_ENVI_FILE,layerstacking,out_name=outname,map_info=map_info,wl=wl

  ENDFOR
  PRINT, '处理完成 : ', SYSTIME()
  ENVI_BATCH_EXIT

END

下图分别是博主处理数据时所使用的数据头文件信息，以及数据类型对照表和对应数据类型的创建函数，hdr用记事本方式打开就可以看头文件信息。