IDL 计算TVDI

介绍请看：http://blog.sina.com.cn/s/blog_764b1e9d0100wdrr.html

源码：

PRO TVDI,NDVI,LST,NBINS,RES
  RES = -1
  
  SZ1 = SIZE(NDVI,/DIMENSIONS)
  SZ2 = SIZE(LST,/DIMENSIONS)
  IF N_ELEMENTS(SZ1) NE 2 OR N_ELEMENTS(SZ2) NE 2 THEN RETURN
  IF TOTAL(ABS(SZ1 - SZ2)) NE 0 THEN RETURN
  IF ~ISA(NDVI,/NUMBER) AND ~ISA(LST,/NUMBER) THEN RETURN
  
  ;FIND DYR WET PTS
  RANGE = [0.0,1.0]
  NPT_FILTER = 200
  PTS = 50
  DRY_PTS = [0]
  WET_PTS = [0]
  
  IDX = FLOOR(NDVI * NBINS)
  FOR I = 0, NBINS DO BEGIN
    PT = WHERE(IDX EQ I,CNT)
    IF(CNT GT NPT_FILTER) THEN BEGIN
      II = SORT(LST[PT])
      DRY_PTS = [ DRY_PTS,PT[II[INDGEN(PTS)]] ]
      WET_PTS = [ WET_PTS,PT[II[-INDGEN(PTS)]] ]
    ENDIF
  END
  
  IF N_ELEMENTS(DRY_PTS) LE 1 THEN RETURN
  
  DRY_PTS = DRY_PTS[1:*]
  WET_PTS = WET_PTS[1:*]
  
  WET = LINFIT(NDVI[WET_PTS],LST[WET_PTS])
  DRY = LINFIT(NDVI[DRY_PTS],LST[DRY_PTS])
  
  GOOD = WHERE(NDVI GE RANGE[0] AND NDVI LE RANGE[1],CNT)
  IF CNT LE 0 THEN RETURN
  
  RES = MAKE_ARRAY(SIZE = SIZE(NDVI))
  RES[GOOD] = (LST - POLY(NDVI[GOOD],WET)) / (POLY(NDVI[GOOD],DRY) - POLY(NDVI[GOOD],WET) )
  
END

说明：

1. 主要是如何求得散点图上下两条边，我的策略是竖着切开散点，分作n个柱子。每个柱子如果总点数大于200，就统计最大最小的50个点作为干湿边的组成部分。

2. sort排序函数，where查找满足条件的元素的下标和个数，linefit线性拟合，poly根据拟合结果和x求得y。具体的多看帮助,那里说的详细,还有例子。