【IDL代码库】利用IDL实现ENVI下计算地形起伏度

在很多具有栅格分析的软件中都没有提供计算地形起伏度的功能，虽然可以根据现有的工具进行组合计算，但是还是比较麻烦。ENVI/IDL的栅格运算功能强大，效率高。本文利用IDL制作了计算地形起伏度的程序，并且能集成在ENVI中使用。

注意事项：将dixingqifudu.pro文件放到ENVI的安装目录\ProgramFiles\ITT\IDLxx\products\enviXX\save_add下，启动ENVI+IDL，出现在ENVI的topographic菜单下：（或者编译成sav文件，只打开ENVI即可。）

源码：

;;===界面居中的函数======================

PRO CENTERTLB, tlb, x, y, NoCenter=nocenter

  COMPILE_OPT StrictArr
 
  geom = WIDGET_INFO(tlb, /Geometry)
 
  IF N_ELEMENTS(x) EQ 0 THEN xc = 0.5 ELSE xc = FLOAT(x[0])
  IF N_ELEMENTS(y) EQ 0 THEN yc = 0.5 ELSE yc = 1.0 - FLOAT(y[0])
  center = 1 - KEYWORD_SET(nocenter)
  ;
  oMonInfo = OBJ_NEW('IDLsysMonitorInfo')
  rects = oMonInfo -> GetRectangles(Exclude_Taskbar=exclude_Taskbar)
  pmi = oMonInfo -> GetPrimaryMonitorIndex()
  OBJ_DESTROY, oMonInfo
 
  screenSize =rects[[2, 3], pmi]
 
  ; Get_Screen_Size()
  IF screenSize[0] GT 2000 THEN screenSize[0] = screenSize[0]/2 ; Dual monitors.
  xCenter = screenSize[0] * xc
  yCenter = screenSize[1] * yc
 
  xHalfSize = geom.Scr_XSize / 2 * center
  yHalfSize = geom.Scr_YSize / 2 * center
 
  XOffset = 0 > (xCenter - xHalfSize) < (screenSize[0] - geom.Scr_Xsize)
  YOffset = 0 > (yCenter - yHalfSize) < (screenSize[1] - geom.Scr_Ysize)
 
  WIDGET_CONTROL, tlb, XOffset=XOffset, YOffset=YOffset
END
;;;================================================

;自动添加envi菜单
PRO dixingqifudu_define_buttons, buttonInfo
  ENVI_DEFINE_MENU_BUTTON, buttonInfo, VALUE = '地形起伏度计算', $\mathrm{uValue=\prime \prime ,𝑢𝑉𝑎𝑙𝑢𝑒=″, event\_pro\; =\text{'}dixingqifudu\text{'},\; \$ REF\_VALUE\; =\; \text{'}Topographic\text{'},position=2end}$

pro qifudu_result, pp, flag, out_name=out_name ;pp是窗口大小，flag是存贮方式
  compile_opt idl2
   envi_batch_init
  pp=pp;窗口大小
  p1=(pp-1)/2;用于建立窗口大小
 envi_select,fid=fid,pos=pos,dims=dims

  if fid eq -1 then begin             

     envi_batch_exit                  

   return                            

 endif                                

  ENVI_FILE_QUERY, fid, bnames=bnames                                       

  nl=dims[4]-dims[3]+1
  ns=dims[2]-dims[1]+1
  nb=size(pos)
  nb=nb[1]
;  print,nb
  tfid=indgen(nb)
 posss=intarr(nb);因为后面生成的都是单波段文件及r_fid，所以都是0表示第一个波段。
  FOR g=0,nb-1 DO BEGIN
    data1 = ENVI_GET_DATA(DIMS=dims, FID=fid , POS=pos[g])
    inherit = envi_set_inheritance(fid, dims, pos[g], /full)
   data2=data1
    FOR i =p1,ns-p1-1 DO BEGIN

     FOR j =p1, nl-p1-1  DO BEGIN

       ;获取一个pp*pp的窗口
       t= data1[[i-p1]:[i+p1],[j-p1]:[j+p1]]
       tmax=max(t,min=tmin)
       data2[i,j]=tmax-tmin
     ENDfor
    ENDFOR
    
     ENVI_WRITE_ENVI_FILE, data2,inherit=inherit ,r_fid=t,out_name=out_name+bnames[pos[g]] ;,/NO_OPEN
;     ENVI_FILE_MNG, id=t, /REMOVE
     tfid[g]=t  
  ENDfor
    if flag eq 1 then begin
      envi_doit, 'cf_doit', fid=tfid, pos=posss, dims=dims, $\mathrm{remove=0,outname=outname,rfid=rfidendifelsebeginenvidoit,\prime cfdoit\prime ,fid=tfid,pos=posss,dims=dims,𝑟𝑒𝑚𝑜𝑣𝑒=0,𝑜𝑢𝑡𝑛𝑎𝑚𝑒=𝑜𝑢𝑡𝑛𝑎𝑚𝑒,𝑟𝑓𝑖𝑑=𝑟𝑓𝑖𝑑𝑒𝑛𝑑𝑖𝑓𝑒𝑙𝑠𝑒𝑏𝑒𝑔𝑖𝑛𝑒𝑛𝑣𝑖𝑑𝑜𝑖𝑡,\prime 𝑐𝑓𝑑𝑜𝑖𝑡\prime ,𝑓𝑖𝑑=𝑡𝑓𝑖𝑑,𝑝𝑜𝑠=𝑝𝑜𝑠𝑠𝑠,𝑑𝑖𝑚𝑠=𝑑𝑖𝑚𝑠, remove=0,\; /IN\_MEMORY,\; r\_fid=r\_fid endelse end}$

;事件响应程序

pro dixingqifudu_event,event
 widget_control,event.top,get_uvalue=pstate

   if TAG_NAMES(event,/STRUCTURE_NAME) eq 'WIDGET_KILL_REQUEST' then begin     
     WIDGET_CONTROL,event.top,/destroy
     ptr_free,pstate
     RETURN
   endif
  
  widget_control,event.id,get_uvalue=uval
  case uval of
   'list':  begin
              (*pstate).chuangkou=fix(event.str)  
              print, fix(event.str)        
;             if event.index eq -1 then a=fix(event.str)
;             help,fix(event.str)
            end
   'file_m':begin
             if event.value eq '内存' then begin
              (*pstate).flag=0
              widget_control,(*pstate).i_base52,map=0
             endif else begin
              (*pstate).flag=1
              widget_control,(*pstate).i_base52,map=1
             endelse
            end 
   'xuanze':begin
             out_name=Dialog_Pickfile(Path=current, /NoConfirm, $
                           Get_Path=path, Filter=['*.*'],title='选择输出文件名')
                 if out_name eq '' then return
                 widget_control,(*pstate).dir_text,set_value=out_name
                 if (*pstate).flag eq 1 then (*pstate).out_name=out_name
            end
   'dir'  :begin             
             widget_control,event.id,get_value=y
             (*pstate).out_name=y
            end
   'queding':begin
   
             if ((*pstate).flag eq 1) and ((*pstate).out_name eq '') then begin
                 tem=DIALOG_MESSAGE('指定文件存储路径')
                 return
             endif
             qifudu_result,(*pstate).chuangkou, (*pstate).flag, out_name=(*pstate).out_name
             envi_batch_exit
            end
   'quxiao':begin
             widget_control,event.top,/DESTROY
             ptr_free,pstate
            end       
    else:
  endcase
end

pro dixingqifudu,event
compile_opt idl2
;  envi,/restore_base_save_files
; envi_select,fid=fid,pos=pos,dims=dims,/band_only
;  print,dims

 

 tlb=widget_base(title='地形起伏度计算',/COLUMN)
  tbase0=widget_base(tlb,/COLUMN,/frame)
  tbase1=widget_base(tbase0,/row)
 lab1=widget_label(tbase1,value='窗口大小（奇数）：   ')
  valu=['3','5','7','9','11']
 list=WIDGET_COMBOBOX(tbase1,value=valu,/EDITABLE,xsize=100,uvalue='list')
 lab2=widget_label(tbase0,value='说明：可键入参数，回车。')
; lab2=widget_label(tbase0,value='',ysize=10)
  lab2=widget_label(tlb,value='',ysize=10)
 
 i_base5=WIDGET_BASE(tlb,/COLUMN,xsize=230,xpad=4,/FRAME)
 k_lable2=widget_label(i_base5,value='选择保存方式：',/ALIGN_LEFT)
  file_m= CW_BGROUP(i_base5, /row, /EXCLUSIVE, /NO_REL, /RETURN_name,UVALUE='file_m',$
    ['文件','内存'], SET_VALUE=0)
 i_base52=WIDGET_BASE(i_base5,xsize=250,/row)
 label3=widget_label(i_base52,value='输出文件',/align_left)
 xuanze=widget_button(i_base52,value='选择',uvalue='xuanze')
 dir_text=WIDGET_TEXT(i_base52,ysize=1,/EDITABLE,UVALUE='dir',/ALL_EVENTS)
  ibase6=widget_base(tlb,/row,/ALIGN_CENTER)
 queding=WIDGET_BUTTON(ibase6,value='确定',uvalue='queding')
 queding=WIDGET_BUTTON(ibase6,value='取消',uvalue='quxiao')
  CENTERTLB, tlb
  widget_control,tlb,/realize
  pstate=ptr_new({dir_text:dir_text,$\mathrm{chuangkou:3,𝑐ℎ𝑢𝑎𝑛𝑔𝑘𝑜𝑢:3, flag:1,out\_name:\text{'}\text{'},\$ i\_base52:i\_base52\},/no\_copy) widget\_control,tlb,set\_uvalue=pstate xmanager,\text{'}dixingqifudu\text{'},tlb,/NO\_BLOCK end}$