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

温度植被干旱指数TVDI(Temperature Vegetation Dryness Index)是一种基于光学与热红外遥感通道数据进行植被覆盖区域表层土壤水分反演的方法。数据来源请引用：地理遥感生态网科学数据注册与出版系统

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