基于 GEE 计算温度植被干旱指数 TVDI

目录

1 温度植被干旱指数（TVDI）

2 完整代码

3 运行结果

---

1 温度植被干旱指数（TVDI）

今天来讲讲在GEE中如何计算温度植被干旱指数（TVDI）。首先，来明确一下有关概念：

“The temperature vegetation dryness index (TVDI) method based on the vegetation index/temperature trapezoid eigenspace (VITT) (Sandholt et al., 2002) also belongs to the category of moisture index method. Thermal inertia method is the approach using thermal infrared remote sensing data to monitor soil moisture.”

译文：“基于植被指数/温度梯形特征空间（VITT）的温度植被干旱指数（TVDI）方法（Sandholt等，2002）也属于湿度指数方法的范畴。热惯性方法是利用热红外遥感数据监测土壤湿度的方法。”

可以发现，TVDI是在NDVI/LST特征空间中构建的，这个特征空间具有一定的分布特征，表现如下图：

NDVI/LST特征空间

NDVI/LST特征空间

NDVI/LST特征空间

而相应地，TDVI的计算公式如下：

式中，Ts为修正后的地表温度；TsMAX为同VI值下最大地表温度，对应VI-Ts干边；TsMIN为同VI值下最小地表温度，对应VI-Ts湿边。

VI和Ts的散点图在植被覆盖和土壤湿度变化较小时为梯形分布，NDVI最高最低温度可通过最大最小值方法进行提取，通过像元进行线性拟合，得出对应的干边、湿边为：

式中，a1、b1为干边拟合方程系数；a2、b2为湿边拟合方程系数；VI为植被干旱指数。

在以往研究中，都是在本地计算TVDI。那么，能不能在GEE中计算TVDI呢？也是可以的，难点在于如何得到干边、湿边的拟合曲线，然后只需要按照TVDI的计算公式就可以计算得到TVDI的全局分布图了。

2 完整代码

var roi = table;
Map.addLayer(roi,{'color':'grey'},'studyArea');
Map.centerObject(roi);

var NDVI = ee.ImageCollection('MODIS/061/MOD13A2')
                  .filter(ee.Filter.date('2020-07-01', '2020-08-01'))
                  .select('NDVI')
                  .mean()
                  .multiply(0.0001)
                  .clip(roi);
var LST = ee.ImageCollection('MODIS/061/MOD11A2')
                  .filter(ee.Filter.date('2020-07-01', '2020-08-01'))
                  .select('LST_Day_1km')
                  .mean()
                  .multiply(0.02)
                  .clip(roi);

var NDVI_LST = NDVI.addBands(LST);
print("NDVI_LST",NDVI_LST);

//采样
var Points = ee.FeatureCollection.randomPoints({          
    region:roi,
    points:3000
    });
// Map.addLayer(Points, {'color':'green'}, 'Points');

var SamplePoint = NDVI_LST.sampleRegions({//提取属性
  collection: Points,
  scale: 1000
});
//print("SamplePoint",SamplePoint);


var SamplePointList = SamplePoint.reduceColumns({
  reducer:ee.Reducer.toList().repeat(2),//通过组合指定数量的给定还原器副本来创建一个还原器。
  selectors:['NDVI','LST_Day_1km']      //输出名称与给定的还原器相同，2列表。
}).get('list');

var SamplePointList2 = ee.List(ee.List(SamplePointList).get(0)).zip(ee.List(SamplePointList).get(1));

var sortColumn = 0;
var keyValues = SamplePointList2.map(function(inner){
    return ee.List(inner).get(sortColumn);
});
//将列表排序为升序。'keyValues'参数被首先排序，'list'中的元素也将以同样的顺序排列。
var sorted = SamplePointList2.sort(keyValues);

var ndviList = sorted.map(function(item){
  return ee.List(item).get(0);
});
var lstList = sorted.map(function(item){
  return ee.List(item).get(1);
});

var tmpList = ee.List.sequence(0,1,0.01);
// print("tmpList",tmpList);
var tmpListSize = tmpList.size().getInfo();

var minList = ee.List([]);
var maxList = ee.List([]);

var tmpList_V2 = tmpList.map(function(item){
  var tmpValue = sorted.map(function(value){
     var listValue = ee.Number(ee.List(value).get(0));
     var listNull = [null,null];
     var upValue = listValue.lt(ee.Number(item).add(0.005));
     var downValue = listValue.gt(ee.Number(item).subtract(0.005));
     var tmpV = ee.Algorithms.If(downValue.and(upValue),
                                value,
                               listNull);
     return tmpV;                             
  });
  return tmpValue;
});


for(var i = 0; i<tmpListSize; i++){
  var currentList = ee.List(tmpList_V2.get(i));
  var currentList_x = currentList.map(function(item){
    return ee.List(item).get(0);
  }).removeAll([null]);
  var currentList_y = currentList.map(function(item){
    return ee.List(item).get(1);
  }).removeAll([null]);
  var sorted_V2 = ee.List(currentList_x).zip(currentList_y).sort(currentList_y);
  var listLength = currentList_y.size().subtract(1);
  // minList = minList.add(sorted_V2.get(0));
  var tmpMin = ee.Algorithms.If(listLength.neq(-1),
                                sorted_V2.get(0),
                                [null,null]);
  var tmpMax = ee.Algorithms.If(listLength.neq(-1),
                                sorted_V2.get(listLength),
                                [null,null]);
  minList = minList.add(tmpMin);
  maxList = maxList.add(tmpMax);
}


var minList_x = minList.map(function(item){
  return ee.List(item).get(0);
}).removeAll([null]);
var minList_y = minList.map(function(item){
  return ee.List(item).get(1);
}).removeAll([null]);

var maxList_x = maxList.map(function(item){
  return ee.List(item).get(0);
}).removeAll([null]);
var maxList_y = maxList.map(function(item){
  return ee.List(item).get(1);
}).removeAll([null]);

var minList = minList_x.zip(minList_y);
print("minList",minList);
var maxList = maxList_x.zip(maxList_y);
print("maxList",maxList);

var linearFit_min = ee.Dictionary(minList.reduce(ee.Reducer.linearFit()));
print("linearFit_min",linearFit_min);
var linearFit_max = ee.Dictionary(maxList.reduce(ee.Reducer.linearFit()));
print("linearFit_max",linearFit_max);

var min_y = ee.Array(ee.List(ndviList).map(function(x) {
  var y = ee.Number(x).multiply(linearFit_min.get('scale')).add(linearFit_min.get('offset'));
  return y;
}));
var max_y = ee.Array(ee.List(ndviList).map(function(x) {
  var y = ee.Number(x).multiply(linearFit_max.get('scale')).add(linearFit_max.get('offset'));
  return y;
}));

var yArr = ee.Array.cat([lstList, min_y,max_y], 1);
print('yArr',yArr);

print(ui.Chart.array.values({
  array: yArr, // lstList yArr
  axis: 0,
  xLabels: ndviList})
  .setChartType('ScatterChart')
  .setOptions({
    legend: {position: 'none'},
    hAxis: {'title': 'NDVI'},
    vAxis: {'title': 'LST'},
    series: {
      0: {
        pointSize: 1,
        // dataOpacity: 0.5,
      },
      1: {
        pointSize: 0,
        lineWidth: 2,
      },
      2: {
        pointSize: 0,
        lineWidth: 2,
      }
    }
  })
);

// apply to Image
var a = linearFit_max.getNumber('offset');
var b = linearFit_max.getNumber('scale');
var c = linearFit_min.getNumber('offset');
var d = linearFit_min.getNumber('scale');
var lst_max = ee.Image(NDVI.multiply(b).add(a));
var lst_min = ee.Image(NDVI.multiply(d).add(c));
var tvdi = (LST.subtract(lst_min)).divide(lst_max.subtract(lst_min));
print(tvdi);

3 运行结果

运行结果

研究区温度植被干旱指数（TVDI）的分布图