TVDI

2021/10/31

# gma.rsvi.TVDI(VI, LST, ValidVIRange = [0.2, 0.8], ValidLSTRange = [200, 350], Delta = 0.01, VINoData = None, LSTNoData = None, OutNoData = -1) 1.0.2 +


功能:【TVDI】。计算温度植被干旱指数。

参数:

 VI:array 。 植被指数。

 LST:array 。 陆表温度(K)。

可选参数:

  ValidVIRange = list 。参与干湿边拟合植被指数的有效值范围,默认 0.2 ~ 0.8。

  ValidLSTRange = list 。参与干湿边拟合陆表温度的有效值范围,默认 200 ~ 350(K)。

  Delta = float 。干湿边拟合过程中植被指数采样的样本间距。默认为 0.01。

  VINoData = float 。植被指数数据中的无效值。默认不设置(None)。

  LSTNoData = float 。陆表温度数据中的无效值。默认不设置(None)。

  OutNoData = float 。TVDI结果中的无效值。默认值为 -1。

注意

VINoData, LSTNoData 的区域在输出数据中都将被改为 OutNoData。

**返回:**类型:namedtuple。包含 TVDI(array),参与拟合的VI(array),湿边LST(array), 干边LST(array),湿边方程系数(flaot),(flaot),干边方程系数(flaot),常数(flaot)。


示例:

from gma import rsvi
1

基于数组类(1 维或多维)

NDVI = [0.5251, 0.5092, 0.4618, 0.4304, 0.4494, 0.4544, 0.4982, 0.6308, 0.5271, 0.4489]
LST = [302.72, 302.98, 303.64, 304.68, 303.7 , 302.94, 302.78, 300.64, 301.98, 302.12]
rsvi.TVDI(NDVI, LST)
1
2
3

>>> TVDI(TVDI=array([ 1.93649901, 1.79681104, 1.2003987 , 2.40119038, 0.79481429, -1.04221242, 0.37934697, -2.2399279 , -1.1514285 , -3.38829638]),
  XVI=array([0.435, 0.445, 0.455, 0.465, 0.495, 0.505, 0.525, 0.635]),
  WetLST=array([304.68, 302.12, 302.94, 303.64, 302.78, 302.98, 301.98, 300.64]),
  DryLST=array([304.68, 303.7 , 302.94, 303.64, 302.78, 302.98, 302.72, 300.64]),
  WetA=-14.890410958904331,
  WetB=310.0907534246576,
  DryA=-16.835616438356006,
  DryB=311.3436301369862)

警告

计算 TVDI 时值过少会引起拟合结果较差或异常!要保证 TVDI 的结果质量需保证足够数量的有效数据!

基于栅格(MODIS VI/LST)

from gma import io
# 读取栅格文件至数据集
NDVISet = io.ReadRaster("MOD13Q1_LY_NDVI_20220407.tif")
LSTSet = io.ReadRaster("MOD11A2_LY_LST_20220407.tif")

# 提取数据集的仿射变换和坐标系
Proj = NDVISet.Projection
Geot = NDVISet.GeoTransform
# 记录不参与拟合和计算的 VI 值和 LST 值( ValidVIRange,ValidLSTRange 范围外的数据不参与拟合)
VINoData = NDVISet.NoData
LSTNoData = LSTSet.NoData

# 依据数据单位进行单位换算
NDVI = NDVISet.ToArray() * 1e-4 
LST = LSTSet.ToArray() * 0.02

TVDI = rsvi.TVDI(NDVI, LST, VINoData = VINoData, LSTNoData = LSTNoData, OutNoData = VINoData)

# 将结果保存为 GTiff 格式
io.SaveArrayAsRaster(TVDI.TVDI, 
                     'MODIS_LY_TVDI_20220407.tif', 
                     Projection = Proj,
                     Transform = Geot,
                     DataType = 'Float32',
                     NoData = VINoData)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

NDVI、LST 和 TVDI 计算结果:

TVDI 干湿边拟合:

import matplotlib.pyplot as plt
from gma import math
# 配置参数
PAR = {'font.sans-serif': 'Times New Roman',
       'axes.unicode_minus': False,
       'mathtext.default': 'default',
      }
plt.rcParams.update(PAR)
plt.figure(figsize = (8, 5), dpi = 300)

# 干边散点图、拟合图和方程
plt.scatter(TVDI[1], TVDI[3], c = 'pink', s = 10)
DryP = TVDI[5]
PolyfitD = TVDI[1] * DryP[0] + DryP[1]
plt.plot(TVDI[1], PolyfitD, c = 'red')
DryR2 = math.Evaluation(PolyfitD, TVDI[3]).R2()
DryInfo = f'Dry: y = {"{:.4f}".format(DryP[0])} * VI + {"{:.4f}".format(DryP[1])}\n{" " * 8}$R^2$ = {DryR2}'
plt.text(0.58, 305, DryInfo, c = 'red', linespacing = 1.5)

# 湿边散点图、拟合图和方程
plt.scatter(TVDI[1], TVDI[2], c = 'lightblue', s = 10)
WetP = TVDI[4]
PolyfitW = TVDI[1] * WetP[0] + WetP[1]
plt.plot(TVDI[1], PolyfitW, c = 'blue')
WetR2 = math.Evaluation.Evaluation(PolyfitW, TVDI[2]).R2()
WetInfo = f'Wet: y = {"{:.4f}".format(WetP[0])} * VI + {"{:.4f}".format(WetP[1])}\n{" " * 8}$R^2$ = {WetR2}'
plt.text(0.58, 291, WetInfo, c = 'blue', linespacing = 1.5)

# 其他配置
plt.grid(True, linestyle = (0,(6,6)), linewidth = 0.4)
plt.ylim(280, 310)
plt.xlabel('VI', loc = 'right', fontsize = 12)
plt.ylabel('LST(K)', loc = 'top', fontsize = 12)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34