加载中…
正文 字体大小:

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

(2018-05-22 10:51:15)
标签:

envi/idl二次开发

tvdi

vtci

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

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

  17.     task=envitask('tvditask')
  18.     task.raster1=raster1
  19.     task.minimum=minimum
  20.     task.step=step
  21.     task.resampling=resampling
  22.     task.display=display
  23.     r=e.ui.selecttaskparameters(task)

  24.     if r eq 'OK' then task.execute else return
  25.   endif
  26.  
  27.   if (TVDI eq 'TVDI') then TVDI=1 else TVDI=0
  28.  
  29.   ;初始化进度条
  30.   Channel = e.GetBroadcastChannel()
  31.   Abort = ENVIAbortable()
  32.   Start = ENVIStartMessage('TVDI-VTCI Task', Abort)
  33.   Channel.Broadcast, Start
  34.   Progress = ENVIProgressMessage('Executing...' , $
  35.     0, Abort)
  36.  
  37.   ;获取两景影像公共区域,并使得结果行列数完全一致
  38.   overLay=GetRasterOverlay  $
  39.     raster1,                  $
  40.     raster2,                  $
  41.     isRaster1Ref=isRaster1Ref,$
  42.     resampling=resampling)

  43.   if n_elements(overLay) ne 2 then return

  44.   if isRaster1Ref then refRaster=overLay[0] else refRaster=overLay[1]
  45.  
  46.   ;初始化一个浮点型Raster存储结果
  47.   rasternew=enviraster     $
  48.     uri=output_Raster,       $
  49.     INHERITS_FROM=refRaster, $
  50.     INTERLEAVE='BSQ'       $
  51.     DATA_TYPE=4)

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

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

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

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

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

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

  76.       return
  77.       BREAK
  78.     ENDIF

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

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

  83.   endfor

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

  87.   NanIndex1=finite(x)
  88.   NanIndex2=finite(y)

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

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

  91.   x=x[xindex]
  92.   y=y[xindex]

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

  97.   _mina=_Min[0]
  98.   _minb=_Min[1]

  99.   _maxa=_Max[0]
  100.   _maxb=_Max[1]

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

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

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

  112.       return
  113.       BREAK
  114.     ENDIF

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

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

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

  122.     rasternew.settile,r,xTiles

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

  131.   rasternew.save
  132.  
  133.   ;进度条结束
  134.   Finish = ENVIFinishMessage(Abort)
  135.   Channel.Broadcast, Finish
  136.  
  137.   ;是否在ENVI中显示结果
  138.   if display then begin
  139.     View = e.GetView()
  140.     Layer = View.CreateLayer(rasternew)
  141.     View.Zoom, /FULL_EXTENT
  142.   endif
  143.  
  144.   ;绘制散点图和报表
  145.   xArr=h['xArr']
  146.   minYarr=h['minYarr']
  147.   maxYarr=h['maxYarr']
  148.   plotWetDry,xArr,minYarr,maxYarr

  149. end

0

阅读 评论 收藏 转载 喜欢 打印举报
已投稿到:
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 不良信息反馈 电话:4006900000 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有