最近写了几个IDL波段分解与合成的代码，分享一下。
一、波段的分解：
1 pro band_decompose 2 compile_opt IDL2 3 4 envi,/Restore_Base_Save_Files 5 envi_batch_init,log_file='batch.txt' 6 7 inputfile='E:L5_12332_20040908 m' 8 envi_open_file,inputfile,r_fid=fid 9 if(fid eq -1) then begin 10 envi_batch_exit 11 return 12 endif 13 14 envi_file_query,fid,dims=dims,nl=nl,ns=ns,data_type=data_type 15 16 ;get the map projection information from the resource file 17 mapinfo=envi_get_map_info(fid=fid) 18 19 ;output to memory and to other format. 20 ;data=envi_get_data(fid=fid,dims=dims,pos=[0]) 21 ;envi_enter_data,data,r_fid=fid 22 ;out_name='E:L5_12332_20040908Bandand1_otherformat.tif' 23 ;envi_output_to_external_format,dims=dims,fid=fid,out_name=out_name,pos=[0],/tiff 24 25 26 27 ;output to ENVI format 28 data=envi_get_data(fid=fid,dims=dims,pos=[7]) 29 30 outputfile='E:L5_12332_20040908Bandand2_1' 31 ;outputfile=data 32 33 OpenW, unit, outputfile, /Get_LUN 34 WriteU, unit, data 35 Free_LUN, unit 36 37 ;create map projection information 38 ;ps=[30,30] 39 ;mc=[342198.231,4489929.046] 40 ;proj=envi_proj_create(/geographic) 41 ;mapinfo=envi_map_info_create($ 42 ;name='Geographic',$ 43 ;mc=mc,$ 44 ;ps=ps,$ 45 ;proj=proj,$ 46 ;/GEOGRAPHIC) 47 48 49 50 ;create the header of file 51 ENVI_SETUP_HEAD, fname=outputfile, $ 52 ns=ns, nl=nl, nb=1, $ 53 interleave=0, data_type=data_type, $ 54 map_info=mapinfo,$ 55 offset=0, /write, /open 56 57 end
可以保持为ENVI格式或者其他格式，代码中有注释。
波段分解批处理：
1 pro band_decompose_batch 2 compile_opt idl2 3 4 envi,/Resotre_Base_Save_Files 5 envi_batch_init,log_file='batch.txt' 6 7 inputfilefolder=dialog_pickfile(/directory,title='Choose the input directory!') 8 ;inputfile=dialog_pickfile(path=EXPAND_PATH('<IDL_DIR>')+PATH_SEP()+'examplesdata',title='打开文件') 9 outputfilefolder=dialog_pickfile(/directory,title='Choose the output directory!') 10 11 inputfilefolderfiles=file_search(inputfilefolder,'*.img',count=num) 12 for j=0,num-1 do begin 13 inputfile=inputfilefolderfiles[j] 14 15 envi_open_file,inputfile,r_fid=fid 16 if(fid eq -1) then begin 17 envi_batch_exit 18 return 19 endif 20 21 envi_file_query,fid,bnames=bnames,dims=dims,nl=nl,ns=ns,nb=nb,data_type=data_type 22 23 mapinfo=envi_get_map_info(fid=fid) 24 25 for i=0,nb-1 do begin 26 data=envi_get_data(fid=fid,dims=dims,pos=[i]) 27 28 ;outputfile=outputfile1+'band_'+strmid(string(i+1),strlen(string(i+1))-1,1) 29 outputfilename1=file_basename(inputfile) 30 outputfilename2=strmid(outputfilename1,0,strlen(outputfilename1)-4) 31 outputfile=outputfilefolder+outputfilename2+'_'+bnames[i] 32 33 OpenW, unit, outputfile, /Get_LUN 34 WriteU, unit, data 35 Free_LUN, unit 36 37 ENVI_SETUP_HEAD, fname=outputfile, $ 38 ns=ns, nl=nl, nb=1, $ 39 interleave=0, data_type=data_type, $ 40 map_info=mapinfo,$ 41 offset=0, /write, /open 42 43 endfor 44 45 endfor 46 47 end
二、波段合成
IDL帮助文档里的一种方法：
1 PRO Band_ENVI_LAYER_STACKING_DOIT 2 3 compile_opt IDL2 4 5 ; 6 7 ; First restore all the base save files. 8 9 ; 10 11 envi, /restore_base_save_files 12 13 ; 14 15 ; Initialize ENVI Classic and send all errors 16 17 ; and warnings to the file batch.txt 18 19 ; 20 21 envi_batch_init, log_file='batch.txt' 22 23 ; 24 25 ; Open the first input file. Also open the 26 27 ; single-band DEM file to layer stack with 28 29 ; this file. 30 31 ; 32 33 envi_open_file, 'E:L5_12332_20040908Band m_band1', r_fid=t_fid 34 35 if (t_fid eq -1) then begin 36 37 envi_batch_exit 38 39 return 40 41 endif 42 43 ; 44 45 ; Open the second input file. 46 47 ; 48 49 envi_open_file, 'E:L5_12332_20040908Band m_band2', r_fid=d_fid 50 51 if (d_fid eq -1) then begin 52 53 envi_batch_exit 54 55 return 56 57 endif 58 59 ; 60 61 ; Use all the bands from both files 62 63 ; and all spatial pixels. First build the 64 65 ; array of FID, POS and DIMS for both 66 67 ; files. 68 69 ; 70 71 envi_file_query, t_fid, $ 72 73 ns=t_ns, nl=t_nl, nb=t_nb 74 75 envi_file_query, d_fid, $ 76 77 ns=d_ns, nl=d_nl, nb=d_nb 78 79 nb = t_nb + d_nb 80 81 fid = lonarr(nb) 82 83 pos = lonarr(nb) 84 85 dims = lonarr(5,nb) 86 87 for i=0L,t_nb-1 do begin 88 89 fid[i] = t_fid 90 91 pos[i] = i 92 93 dims[0,i] = [-1,0,t_ns-1,0,t_nl-1] 94 95 endfor 96 97 ; 98 99 for i=t_nb,nb-1 do begin 100 101 fid[i] = d_fid 102 103 pos[i] = i-t_nb 104 105 dims[0,i] = [-1,0,d_ns-1,0,d_nl-1] 106 107 endfor 108 109 ; 110 111 ; Set the output projection and 112 113 ; pixel size from the TM file. Save 114 115 ; the result to disk and use floating 116 117 ; point output data. 118 119 ; 120 121 out_proj = envi_get_projection(fid=t_fid, $ 122 123 pixel_size=out_ps) 124 125 out_name = 'C:Band1_2' 126 127 out_dt = 4 128 129 ; 130 131 ; Call the layer stacking routine. Do not 132 133 ; set the exclusive keyword allow for an 134 135 ; inclusive result. Use cubic convolution 136 137 ; for the interpolation method. 138 139 ; 140 141 envi_doit, 'envi_layer_stacking_doit', $ 142 143 fid=fid, pos=pos, dims=dims, $ 144 145 out_dt=out_dt, $ 146 147 out_name=out_name, $ 148 149 interp=2, $ 150 151 out_ps=out_ps, $ 152 153 out_proj=out_proj, r_fid=r_fid 154 155 ; 156 157 ; Exit ENVI Classic 158 159 ; 160 161 ;envi_batch_exit 162 163 END
另一种方法：
1 pro band_envi_layer_stacking_doit_test 2 compile_opt idl2 3 4 envi,/restore_base_save_files 5 envi_batch_init,log_file='batch.txt' 6 7 inputfiles=['E:L5_12332_20040908Band m_band1','E:L5_12332_20040908Band m_band2'] 8 outputfile='c: m_band_12' 9 10 fids=lonarr(n_elements(inputfiles)) 11 dimses=lonarr(5,n_elements(inputfiles)) 12 poses=lonarr(n_elements(inputfiles)) 13 14 for i=0,n_elements(inputfiles)-1 do begin 15 envi_open_file,inputfiles[i],r_fid=fids1 16 envi_file_query,fids1,ns=ns,nl=nl,nb=nb 17 18 fids[i]=fids1 19 dimses[0,i]=[-1,0,ns-1,0,nl-1] 20 proj=envi_get_projection(fid=fids,pixel_size=out_ps) 21 22 endfor 23 24 envi_doit,'envi_layer_stacking_doit',$ 25 fid=fids,pos=poses,dims=dimses,$ 26 out_dt=4,out_name=outputfile,$ 27 interp=2,out_ps=out_ps,$ 28 out_proj=proj,r_fid=r_fid 29 30 end