基于MODIS数据的NDVI与LST相关性分析（IDL代码实现）

最新推荐文章于 2024-09-30 19:00:07 发布

遇见飖雪

最新推荐文章于 2024-09-30 19:00:07 发布

阅读量7.4k

点赞数 11

分类专栏： IDL 可视化开发语言遥感

本文链接：https://blog.csdn.net/qq_41799597/article/details/106284888

版权

遥感同时被 3 个专栏收录

6 篇文章 0 订阅

订阅专栏

IDL

5 篇文章 1 订阅

订阅专栏

可视化开发语言

5 篇文章 0 订阅

订阅专栏

1 数据预处理

（1）数据提取
我们可以选取2018年5月初华东地区MODIS中的MOD11A2和MOD13A2的16天合成LST和NDVI产品数据，下载地址：MODIS数据下载
网站下载数据需要注册，在此声明，下载数据不需要“科学上网”，但注册需要谷歌Gmail邮箱，注册Gmail邮箱需要“科学上网”。以下是我们下载数据，数据下载中可以通过Python脚本的方法对数据进行批量下载，由于我们只需要少量数据，不必批量下载，所以在此不加赘述。以下是我们的原始数据在ENVI中的显示：
图1 原始数据 Modis的NDVI产品数据

图2 原始数据 Modis的LST产品数据

由于MODIS产品数据下载之后不能直接使用，因此需要对数据进行预处理。

首先EarthData网站上获取的数据为*.hdf格式，此时不能直接有ENVI软件直接处理，需要批量提取，这时就要用到NASA提供的MODIS Reprojection Tool，此工具虽不能实现全自动的批量提取，但是可以实现按月进行数据的提取及拼接，以下是利用MRT软件提取数据参数设置：
图3 MRT软件提取MODIS产品参数设置
NDVI产品数据也用以上同样的方法输入MOD11A2数据提取，输出LST和NDVI数据。

（2）数据复原
此时数据依然不能直接使用，由于LST数据是16bit的，其显示范围是0－65535，有效范围则是7500－65535，需要乘以0.02才能得到开尔文温度。而NDVI数据为了减少储存空间，也是正常值的10000倍，所以应该给其乘以0.0001得到正常值，以上都是利用EMVI的Bandmath工具进行，在此不加赘述。
图4 数据复原统计NDVI

l图5 数据复原统计LST

（3）数据裁剪
由于涉及范围比较广，所以我们对数据进行裁剪处理，裁取山东地区的一小部分作为分析样本进行数据回归分析与相关性分析，具体范围如下专题图：
图6 裁剪后的NDVI产品彩色分级显示
图7 裁剪后的LST产品彩色分级显示

2 代码与步骤

Pro Regress_ndvi_lst

;input the image

fn_ndvi=**dialog_pickfile**(title='Please Input the NDVI Image:')

fn_lst=**dialog_pickfile**(title='Please Input the LST Image:')

;read the image

ndvi_img=**read_image**(fn_ndvi)

lst_img=**read_image**(fn_lst)

;get the size of image

sz_ndvi=**size**(ndvi_img)

;get the image's numbers of columns and rows

ndvi_columns=sz_ndvi[**1**] & ndvi_rows=sz_ndvi[**2**]

;get the size of image

sz_lst=**size**(lst_img)

;get the image's numbers of columns and rows

lst_columns=sz_lst[**1**] & lst_rows=sz_lst[**2**]

;samples

random_c=**randomu**(seed,ndvi_columns)

random_r=**randomu**(seed,ndvi_rows)

sort_c=**sort**(random_c)

sort_r=**sort**(random_r)

ndvi_samples=ndvi_img[sort_c[**0**:**10**],*]

lst_samples=lst_img[sort_c[**0**:**10**],*]

ndvi_get=**reform**(ndvi_samples,**11***ndvi_rows, **1**)

lst_get=**reform**(lst_samples,**11***ndvi_rows, **1**)

ref_ndvi=**reform**(ndvi_get)

ref_lst=**reform**(lst_get)

;ref_ndvi1=reform(ndvi_img,ndvi_columns*ndvi_rows,1)

;ref_lst1=reform(lst_img,lst_columns*lst_rows,1)

;ref_ndvi=ref_ndvi1[0:100]

;ref_lst=ref_lst1[0:100]

;regress analysis

fit_ndvi_lst=**regress**(ref_ndvi,ref_lst,const=b,correlation=r,yfit=lst_estimated)

;plot points picture

p1=**plot**(ref_ndvi,ref_lst,xtitle='NDVI',ytitle='LST(K)',dimensions=[**600**,**400**], symbol=**1**,color='red',sym_size=**1.0**,linestyle=**6**,$

   title='The Regress of NDVI & LST',window_title='NDVI & LST',$

    xrange=[**0**,**1**],yrange=[**290**,**310**])

;set the range of x and y axis

x=ref_ndvi

y=fit_ndvi_lst[**0**]*x+b

fit=**plot**(x,y,linestyle=**1**,color='blue',thick=**2**,/current,/overplot)

p2=**plot**([**0**,**1**],[**290**,**310**],linestyle=**6**,/current,/overplot)

str_equation='Y='+**string**(fit_ndvi_lst[**0**],format='(f6.2)')+'*X+'+**string**(b,format='(f6.2)')

str_correlation='R='+**string**(r,format='(f5.2)')

t1=**text**(**0.20**,**0.80**,str_equation,font_size=**12**,target=p1)

t2=**text**(**0.20**,**0.75**,str_correlation,font_size=**12**,target=p1)

;o_fn='Scatter_ndvi_lst.emf'

;p1.save,o_fn,border=40

;p2.save,o_fn,border=40

;**********

;plot points picture

p3=**plot**(ref_lst,lst_estimated,xtitle='LST(K)',ytitle='LST-Estimated(K)',dimensions=[**600**,**400**],$

    symbol=**24**,color='red',sym_size=**1.0**,title='LST & LST-Estimated',$

    window_title='Cal_LST_Linear',linestyle=**6**,xrange=[**290**,**310**],yrange=[**290**,**310**])

;set the range of x and y axis

p4=**plot**([**290**,**310**],[**290**,**310**],thick=**2**,color='blue',/current,/overplot)

MAE=**mean**(**abs**(ref_lst-lst_estimated))

RMSE=**sqrt**(**mean**((ref_lst-lst_estimated)^**2**))

str_MAE='MAE='+**string**(MAE,format='(f5.2)')

str_RMSE='RMSE='+**string**(RMSE,format='(f5.2)')

t3=**text**(**0.20**,**0.80**,str_MAE,font_size=**12**,target=p3)

t4=**text**(**0.20**,**0.75**,str_RMSE,font_size=**12**,target=p3)

;o_fn='Scatter_lst_lst_estimated.emf'

;p3.save,o_fn,border=40

;p4.save,o_fn,border=40

**End**

3 结果与分析：

由于代码中是随机数采样，每次采样训练样本数据不同，所以每次运行代码可以得到不同的模型。以下是模型显示：
模型1:
图8 模型1线性回归

图9 模型1精确度与误差分析

模型2:
图11 模型2精确度与误差分析
图10 模型2线性回归

模型3:
图13 模型3精确度与误差分析

图12 模型3线性回归

也可以修改代码中的采样数目，来改变训练样本数目，模型也会发生变化。如将其中的[0:10]改为[0:20],对应的11修改为21就可增加训练样本数目：

代码：

ndvi_samples=ndvi_img[sort_c[**0**:**20**],*]

lst_samples=lst_img[sort_c[**0**:**20**],*]

ndvi_get=**reform**(ndvi_samples,**21***ndvi_rows, **1**)

lst_get=**reform**(lst_samples,**21***ndvi_rows, **1**)

模型4:

图14 模型4线性回归

图15 模型3精确度与误差分析

结论分析：

相关系数（R），是衡量两个变量之间相关程度的系数，是判定变量之间线性相关性的一个相对指标。相关系数用字母R表示，相关系数R取值在±1之间，当R为0时，表示两个变量绝对不相关；当R大于0时，两个变量正相关，即你增加我也增加，你减少我也减少；当R小于0时，两个变量负相关，即你增加我减少，你减少我增加；当R等于1或-1时，表示两个变量绝对相关。

相关系数R越接近于±1，两个变量之间相关性越强。一般认为：当R的绝对值为0.7或更大时，两个变量高度相关，即强相关；当R的绝对值在0.5-0.7之间时，两个变量中度相关；当R绝对值在0.3-0.5之间时，两个变量弱相关；当R绝对值低于0.3时，说明两个变量之间几乎不存在相关关系。

平均绝对误差（MAE），MAE 的值越小，说明预测模型拥有更好的精确度。范围[0,+∞)，当预测值与真实值完全吻合时等于0，即完美模型；误差越大，该值越大。

均方根误差（RMSE），它表示误差的平方的期望值，范围[0,+∞)，当预测值与真实值完全吻合时等于0，即完美模型；误差越大，该值越大。

综上，以上四种模型的相关系数R的绝对值值均在0.5-0.7之间，说同一时间NDVI值与LST中度相关。其中第四种模型R的绝对值最大，相关性更高，进一步说明，采样数据越多模型相关性越好；对比四种模型的平均绝对误差（MAE）和均方根误差（RMSE），第四种模型这两个量较前三种较小，说明其拥有更好的精确度，且误差较小，更进一步说明，样本数据越多，其拟合度越好，模型越精确。