【IDL代码库】一个完整的ENVI扩展工具源码

以TVDI VTCI扩展工具为例,为广大遥感爱好者提供一个完整ENVI/IDL二次开发示例,包括算法编写、数据分块处理、绘图、IDL界面搭建、事件响应和自定义ENVITask等内容。
完整源码下载:https://pan.baidu.com/s/1REgEB0R3yWOeLIRASi44xg 提取码:cuk6
以下仅贴出关键代码,作者水平有限,难免有不足,还望指正。
;+
; :Author: Hanzt
; :Description
; 通过拟合干、湿边方程,求TVDI或VTCI指数
;-
;关键代码
pro TVDItask,                $
 raster1=raster1,           $
 raster2=raster2,           $
 TVDI=TVDI,                 $
 minimum=minimum,           $
 isRaster1Ref=isRaster1Ref,  $
 step=step,                 $
 resampling=resampling,     $
 output_Raster=output_Raster,$
 display=display

  COMPILE_OPT idl2
  e =envi(/current)
 
  if raster1.uri eq raster2.uri then begin
   void=DIALOG_MESSAGE('The input NDVI and LST raster must be different!',$
     /info)

   task=envitask('tvditask')
   task.raster1=raster1
   task.minimum=minimum
   task.step=step
   task.resampling=resampling
   task.display=display
   r=e.ui.selecttaskparameters(task)

   if r eq 'OK' then task.execute else return
  endif
 
  if (TVDI eq 'TVDI') then TVDI=1 else TVDI=0
 
  ;初始化进度条
  Channel = e.GetBroadcastChannel()
  Abort = ENVIAbortable()
  Start = ENVIStartMessage('TVDI-VTCI Task', Abort)
  Channel.Broadcast, Start
  Progress = ENVIProgressMessage('Executing...' , $
   0, Abort)
 
  ;获取两景影像公共区域,并使得结果行列数完全一致
  overLay=GetRasterOverlay(  $
   raster1,                 $
   raster2,                 $
   isRaster1Ref=isRaster1Ref,$
   resampling=resampling)

  if n_elements(overLay) ne 2 then return

  if isRaster1Ref then refRaster=overLay[0] else refRaster=overLay[1]
 
  ;初始化一个浮点型Raster存储结果
  rasternew=enviraster(     $
   uri=output_Raster,      $
   INHERITS_FROM=refRaster, $
   INTERLEAVE='BSQ',       $
   DATA_TYPE=4)

  IF overLay[0].metadata.HasTag('data ignore value') THEN $
   dataIgnoreValue=overLay[0].metadata['data ignore value']

  ;设置分块大小
  tileSize=[2048,2048]
 
  tileSize[0] = overLay[0].ncolumns gt tileSize[0] ? tileSize[0] : overLay[0].ncolumns
  tileSize[1] = overLay[0].nRows gt tileSize[1] ? tileSize[1] : overLay[0].nRows

  xTiles=overLay[0].CreateTileIterator(tile_Size=tileSize)
  yTiles =overLay[1].CreateTileIterator(tile_Size=tileSize)

  x=!null
  y=!null
 
  ;两次分块
  ;第一次获取参与拟合的NDVI和LST数组
  ;第二次分块计算TVDI或VTCI指数
 
  ;获取参与拟合的NDVI和LST数组
  for i=0,xTiles.ntiles-1 do begin

   percentProgress = Round(i* 100.0/(xTiles.ntile*2))
   Progress.Percent = percentProgress
   Channel.Broadcast, Progress

   IF (Abort.Abort_Requested) THEN BEGIN
     Finish = ENVIFinishMessage(Abort)
     Channel.Broadcast, Finish

     return
     BREAK
   ENDIF

   xTmp=xTiles.next()
   yTmp=yTiles.next()

   x=[x,xTmp[0:*:step]]
   y=[y,yTmp[0:*:step]]

  endfor

  if isa(dataIgnoreValue,/number) then $
   xIndex=where(x ne dataignorevalue and x gt minimum) else $
   xIndex=where( x gt minimum)

  NanIndex1=finite(x)
  NanIndex2=finite(y)

  nanIndex=where(NanIndex1 eq 0 or NanIndex2 eq 0)

  if nanIndex[0] ne -1 then xindex=setdifference(xIndex,nanIndex)

  x=x[xindex]
  y=y[xindex]

  ;获取回归参数
  h=para_tvdi(x,y,minimum=minimum)
  _Min=h['Min']
  _Max=h['Max']

  _mina=_Min[0]
  _minb=_Min[1]

  _maxa=_Max[0]
  _maxb=_Max[1]

  xTiles.reset
  yTiles.reset
 
  ;求得回归参数后分块计算TVDI或VTCI
  for i=0,xTiles.ntiles-1 do begin

   percentProgress = Round((i+xtiles.ntiles)* 100.0/(xTiles.ntile*2))
   Progress.Percent = percentProgress
   Channel.Broadcast, Progress

   IF (Abort.Abort_Requested) THEN BEGIN
     Finish = ENVIFinishMessage(Abort)
     Channel.Broadcast, Finish

     return
     BREAK
   ENDIF

   x=xTiles.next()
   y=yTiles.next()

    r=cal_TVDI_VTCI(_mina,_minb,_maxa,_maxb,x,y,TVDI=TVDI)

   if isa(dataIgnoreValue,/NUMBER) then begin
     index=where(x eq dataIgnoreValue)
     if index[0] ne -1 then r[index]=dataIgnoreValue
   endif

   rasternew.settile,r,xTiles

  endfor
 
  ;更新结果波段名称
  if TVDI then bnames=['TVDI'] else bnames=['VTCI']
  
  IF rasternew.metadata.HasTag('BAND NAMES') then $
   rasternew.metadata.UpdateItem,'BAND NAMES',bnames ELSE $
   rasternew.metadata.Additem,'BAND NAMES',bnames

  rasternew.save
 
  ;进度条结束
  Finish = ENVIFinishMessage(Abort)
  Channel.Broadcast, Finish
 
  ;是否在ENVI中显示结果
  if display then begin
    View = e.GetView()
    Layer = View.CreateLayer(rasternew)
   View.Zoom, /FULL_EXTENT
  endif
 
  ;绘制散点图和报表
  xArr=h['xArr']
  minYarr=h['minYarr']
  maxYarr=h['maxYarr']
  plotWetDry,xArr,minYarr,maxYarr

end

 

 
posted @ 2022-05-27 17:58  ENVI-IDL技术殿堂  阅读(1495)  评论(0编辑  收藏  举报