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


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