Python计算温度植被干旱指数（TVDI）

温度植被干旱指数TVDI(Temperature Vegetation Dryness Index)是一种基于光学与热红外遥感通道数据进行植被覆盖区域表层土壤水分反演的方法。

1、原理

LST-NDVI特征空间1：

温度植被干旱指数(TVDI)的计算方法为：

TVDI= LST max −LST min LST−LST min

LSTmin=a+b×NDVI

LSTmax=c+d×NDVI

其中
a 、 b 、 c 、 d
为干、湿边拟合系数；

2、Python代码

import gdal
from gdalconst import *
import numpy as np
from glob import glob
from os import path as osp
import os, subprocess
import matplotlib.pyplot as plt# 获取lst、ndvi数据
def get_data(file_ndvi,file_lst):ndvi_tif=gdal.Open(file_ndvi,GA_ReadOnly)lst_tif=gdal.Open(file_lst,GA_ReadOnly)ndvi_band=ndvi_tif.GetRasterBand(1)ndvi=ndvi_band.ReadAsArray()lst_band=lst_tif.GetRasterBand(1)lst=lst_band.ReadAsArray()    return ndvi,lst# 获取投影等信息，用于保存TVDI结果
def get_info(file_ndvi):ndvi_tif=gdal.Open(file_ndvi,GA_ReadOnly)ndvi_band=ndvi_tif.GetRasterBand(1)gt = ndvi_tif.GetGeoTransform()proj = ndvi_tif.GetProjectionRef()dtype = ndvi_band.DataTypereturn gt,proj,dtype# 计算lst的最小值（湿边）和最大值（干边）
def get_min_max(ndvi,lst):MiniList = []MaxList = []# 创建ndvi向量（0到1，间距为0.01）ndvi_vector = np.round(np.arange(0.01, 1.01, 0.01), 2)# 首先找到相同ndvi的lst值for val in ndvi_vector:lst_lst_val = []row, col = np.where((ndvi >= val-0.001) & (ndvi <= val+0.001))# 根据这些ndvi的位置，我们取温度值对应这些位置（行和列）for i in range(len(row)):if np.isfinite(lst[row[i], col[i]]):lst_lst_val += [lst[row[i], col[i]]]# 如果所需的ndvi有lst值，则计算最大值和最小值if lst_lst_val != []:lst_min_val = np.min(lst_lst_val)lst_max_val = np.max(lst_lst_val)else:lst_min_val = np.nanlst_max_val = np.nan# 找到的值被添加到MiniList和MaxList列表中MiniList += [lst_min_val]MaxList  += [lst_max_val]return MiniList,MaxListdef fit(MiniList,MaxList):ndvi_vector = np.round(np.arange(0.01, 1.01, 0.01), 2)MiniList_fin = []ndvi_fin = []for i, val in enumerate(MiniList):if np.isfinite(val):MiniList_fin += [val]ndvi_fin += [ndvi_vector[i]]MinPfit = np.polyfit(ndvi_fin[14:89], MiniList_fin[14:89], 1)MaxList_fin = []ndvi_fin = []for i, val in enumerate(MaxList):if np.isfinite(val):MaxList_fin += [val]ndvi_fin += [ndvi_vector[i]]MaxPfit = np.polyfit(ndvi_fin[14:89], MaxList_fin[14:89], 1)return MinPfit,MaxPfitdef plot_scatter(ndvi,lst,MiniList,MaxList,MinPfit,MaxPfit,scatter_file=None):ndvi_vector = np.round(np.arange(0.01, 1.01, 0.01), 2)a1, b1 = MaxPfita2, b2 = MinPfitlinhamax = [b1 + (a1 * 0), b1 + (a1 * 1)]linhamin = [b2 + (a2 * 0), b2 + (a2 * 1)]plt.plot(ndvi.ravel(), lst.ravel(), "+", color='dimgray', markersize=4)plt.plot(ndvi_vector[14:89], MiniList[14:89], '+', color='b')plt.plot(ndvi_vector[14:89], MaxList[14:89], '+', color='r')if a1>0:plt.plot([0, 1], linhamax, color='r', markersize=8,\label=f"Tsmax = {'%.1f'% b1} + {'%.1f' % abs(a1)} * ndvi")else:plt.plot([0, 1], linhamax, color='r', markersize=8,\label=f"Tsmax = {'%.1f'% b1} - {'%.1f' % abs(a1)} * ndvi")if a2>0:plt.plot([0, 1], linhamin, color='b', markersize=8,\label=f"Tsmin = {'%.1f' % b2} + {'%.1f' % abs(a2)} * ndvi")else:plt.plot([0, 1], linhamin, color='b', markersize=8,\label=f"Tsmin = {'%.1f' % b2} - {'%.1f' % abs(a2)} * ndvi")plt.legend(loc='upper right', fontsize=12)plt.ylim(top=340,bottom=270)plt.xlabel("ndvi")plt.ylabel("lst (K)")plt.title("ndvi vs lst Scatterplot")if scatter_file is not None:plt.savefig(scatter_file)plt.show()def show_tvdi(tvdi,fig_file=None):plt.imshow(tvdi,cmap= 'jet_r',vmax=1,vmin = 0)plt.axis('off')plt.colorbar()plt.title("tvdi")if fig_file is not None:plt.savefig(fig_file)plt.show()def compute_tvdi(ndvi,lst,MinPfit,MaxPfit):a1, b1 = MaxPfita2, b2 = MinPfitTs_max = b1 + (a1 * ndvi)Ts_min = b2 + (a2 * ndvi)TVDI = (lst - Ts_min) / (Ts_max - Ts_min)return TVDIdef save_tvdi(TVDI,gt,proj,dtype,file_out):fname_out   = file_outdriver      = gdal.GetDriverByName('GTiff')data_type   = dtypedset_output = driver.Create(fname_out, TVDI.shape[1], TVDI.shape[0], 1, gdal.GDT_Float32)dset_output.SetGeoTransform(gt)dset_output.SetProjection(proj)dset_output.GetRasterBand(1).WriteArray(TVDI)dset_output.FlushCache()dset_output = Nonedef main(ndvi_file,lst_file,tvdi_file,scatter_file=None,fig_file=None):'''Parameters----------ndvi_file : the file of ndvilst_file : the file of lsttvdi_file : the file use to save tvdiscatter_file : the file use to save scatterfig_file : the file use to save the figure of tvdi'''# 获取ndvi和lst数据ndvi,lst=get_data(ndvi_file,lst_file)ndvi[ndvi<0]=np.nanlst[lst<250]=np.nan# 获取lst的最小值（湿边）和最大值（干边）MiniList,MaxList=get_min_max(ndvi,lst)# 计算tvdi,并保存MinPfit,MaxPfit=fit(MiniList,MaxList)tvdi=compute_tvdi(ndvi,lst,MinPfit,MaxPfit)gt,proj,dtype=get_info(ndvi_file)save_tvdi(tvdi,gt,proj,dtype,tvdi_file)# 显示散点图plot_scatter(ndvi,lst,MiniList,MaxList,MinPfit,MaxPfit,scatter_file)# 展示tvdi结果show_tvdi(tvdi,fig_file)if __name__ == '__main__':ndvi_file=r'*.tif'lst_file=r'*.tif'tvdi_file=r'*.tif'main(ndvi_file,lst_file,tvdi_file)

3、结果展示
LST-NDVI的散点图:

TVDI展示：

如果对您有用的话，别忘了给点个赞哦_ ！数据来源请引用：地理遥感生态网科学数据注册与出版系统

4、参考文献
[1] Du L, Song N, Liu K, et al. Comparison of Two Simulation Methods of the Temperature Vegetation Dryness Index (TVDI) for Drought Monitoring in Semi-Arid Regions of China[J]. Remote Sensing, 2017, 9(2): 177.