• IDL波段分解与合成源代码


    最近写了几个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
  • 相关阅读:
    MongoDB对比关系型数据库
    Swagger 打开时自动折叠
    更改Linux定时任务crontab启动基目录
    linux系统/etc/init.d目录下的开机自启脚本
    vue 中新窗口打开vue页面 (this.$router.resolve)
    树莓派4B如何手动固定IP地址
    树莓派无显示器设置WiFi、开启ssh、开启VNC
    递归
    学习-HTML5
    只是为了表示我还有在敲代码
  • 原文地址:https://www.cnblogs.com/zhzhx/p/3141031.html
Copyright © 2020-2023  润新知