【IDL程式碼庫】一個完整的ENVI擴充套件工具原始碼
阿新 • • 發佈：2022-05-27
以TVDI VTCI擴充套件工具為例，為廣大遙感愛好者提供一個完整ENVI/IDL二次開發示例，包括演算法編寫、資料分塊處理、繪圖、IDL介面搭建、事件響應和自定義ENVITask等內容。
擴充套件工具見：http://blog.sina.com.cn/s/blog_764b1e9d0102y95q.html。
完整原始碼下載：https://pan.baidu.com/s/1REgEB0R3yWOeLIRASi44xg 提取碼：cuk6
以下僅貼出關鍵程式碼，作者水平有限，難免有不足，還望指正。
;+
; :Author: Hanzt
; :Email: [email protected]
; :Description
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