【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 = '地形起伏度计算', $uValue = '',$
    event_pro ='dixingqifudu', $
    REF_VALUE = 'Topographic',position=2
end

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, $remove=0, out_name=out_name, r_fid=r_fid     endif else begin       envi_doit, 'cf_doit', fid=tfid, pos=posss, dims=dims,$
      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, $chuangkou:3,$
   flag:1,out_name:'',$
   i_base52:i_base52},/no_copy)
  widget_control,tlb,set_uvalue=pstate
  xmanager,'dixingqifudu',tlb,/NO_BLOCK

end