实验目的
熟练遥感数据的处理流程；通过探讨地面测量的叶面积指数与遥感观测的植被指数关系，理解地面测量数据与遥感观测数据的联系。
实验内容
- 预处理遥感数据，得到WDRVI指数影像，并提取地面观测LAI对应时间，卫星观测的试验田所在位置的WDRVI均值。
- 对比LAI与WDRVI指数，构建模型，绘制相应的图表，如散点图
- 撰写报告，描述实验过程和结果
实验数据和环境
MODIS卫星观测数据
MODIS数据产品概述包括4级产品：1级产品指L1A数据，已经被赋予定标参数；2级产品：指L1B数据，经过定标定位后的数据，本系统产品是国际标准的EOS-HDF格式；3级产品是在1B数据的基础上，对有遥感器成像过程产生的边缘畸变（Bowtie效应）进行校正，产生L3级产品；4级产品由参数问价提供的参数，对图像进行几何纠正、辐射校正，使图像的每一点都有精确的地理编码、反射率和辐射率。
本次实验所采用的MODIS产品为MOD09Q1，MOD09Q1提供了1-2波段250米8天合成的栅格化3级数据产品。
田间测量的叶面积指数LAI
叶面积指数（Leaf Area Index）是指单位面积土地上所有植物叶片面积的总 和。LAI 对陆地生态系统和大气间能量、水汽和 CO2 交换具有重要影响。通过遥感可以在更大时空尺度上反演 LAI 时空分布。遥感反演 LAI 方法有: (1)植被指数法；(2)像元成分非混合法；(3)直接模式转换 法；(4)间接模式转换法。
本实验中的叶面积指数由实地测量获得，如下图所示。
实验过程
遥感数据预处理
WDRVI与LAI的时间匹配
MOD09Q1产品是间隔八天合成的数据，为了进行较为精确的比较分析，需要根据产品DOY的日期字段与LAI测量时间进行匹配，时间匹配这一环节主要基于python实现，匹配过程遵循以下两个准则：
- 如果LAI观测时间在[DOY,DOY+8-1]这一区间内，则匹配该区间
- 如果不在观测时间内，则选取LAI与DOY日期最近的那一时期影像
提取试验田对应像元的WDRVI值
使用GDAL来提取试验田对应像元的WDRVI值，主要思路是将经纬度也即地理坐标系坐标转为投影坐标系坐标，最后转为图上坐标，进而获取对应像元的WDRVI值，代码见附录。
GDAL的介绍如下：
构建植被指数时间序列
根据所得结果绘制初步的叶面积指数LAI时序变化曲线和WDRVI植被指数时序变化曲线如下所示。研究区玉米地的LAI呈现明显的先增长后减小的总体趋势，且随着不同年份而具有周期性变化规律，非常符合玉米从发芽到枯萎的完整生长情况。同时也可以看出叶面积在不同年份的变化趋势并不是完全一致的，由当年的气候、生长情况等影像而有不同的小型起伏和变化速率。
从下图的两个植被指数时间序列曲线对比，可发现WDRVI时序变化与 LAI 大致相同，但是 WDRVI 在某些时间段呈现多个峰值，主要是因为遥感数据受到云雨等噪声影响，因此需要对时序数据去噪重建时序数据。
原始WDRVI时序数据：
LAI时序数据：
WDRVI时间序列重建
滤波算法能有效消除序列中存在的噪声,在以滤波算法为基础进行NDVI时间序列重建的算法中,基于非对称高斯函数拟合方法、傅立叶拟合系列方法以及S-G滤波拟合方法是3种主要算法。Savitzky-Golay滤波器（通常简称为S-G滤波器）是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变。S-G滤波方法使用滑动窗口对原始植被指数时序数据进行平滑滤波,它在消除噪声、获取长时间序列的变化趋势以及局部突变信息具有较好的效果。
针对提取出来的WRNDVI进行SG滤波一维时间序列重建，根据经验设置滤波参数如下：多项式拟合阶数为3 滑动窗口所包含数据点为7。主要通过调用Python的scipy库中savgol_filter函数来实现，实现代码见附录，平滑结果如图。对比平滑前后的WDRVI时间序列曲线，可以看出SG滤波在保持原有的较多程度的信息的基础上，减小了噪声，最明显的变化体现在不同年份之间的谷值被适当提高。
LAI 与 WDRVI模型建立
基本统计量的计算
使用Seaborn计算并绘制基本统计量的图像，使用pairplot展现变量两两之间的关系（线性或非线性，有无较为明显的相关关系）。其中对角线上为LAI和WDRVI的散点图，其余位置为LAI和WDRVI的数据分布直方图。
相关性分析
为了验证WDRVI和LAI两者之间有较大的关联，首先进行皮尔逊相关分析并绘制散点图，得出两者的相关系数r为0.949，且t=0< 0.05，t检验通过。因此两者存显著的相关性，探究两者的数学表达关系是有依据有意义的。皮尔逊相关系数矩阵如下所示：
线性回归
使用Python的sklearn函数库中的线性回归函数来拟合WDRVI和LAI：
模型的基本参数如下：
| 回归模型参数与R2 | 回归模型参数值与R2值 | ||||
|---|---|---|---|---|---|
| R2 | 0.836 | ||||
| 截距 | -0.844 | ||||
| 斜率 | 0.064 |
最终建立的模型为：
附录
导入函数库
# -*- coding: utf-8 -*-
"""
PROJECT_NAME: RS_Toolbox
FILE_NAME: get_value_by_location
AUTHOR: welt
E_MAIL: tjlwelt@foxmail.com
DATE: 2023/3/28
"""
import os
import os.path as path
import matplotlib.pyplot as plt
import numpy as np
import xlrd
from osgeo import gdal
from osgeo import osr
from scipy.signal import savgol_filter
from sklearn.linear_model import LinearRegression
from scipy.stats import stats
import seaborn as snsimport pandas as pd
获取指定位置的值
def get_img_value(wdrvi_path, lon_lat):
"""
:param wdrvi_path: 图像地址
:param lon_lat: 要获取位置的经纬度
:return: 获取的值
"""
# 读入数据
dataset = gdal.Open(wdrvi_path)
img = dataset.ReadAsArray()
pcs = osr.SpatialReference()
pcs.ImportFromWkt(dataset.GetProjection())
gcs = pcs.CloneGeogCS()
extend = dataset.GetGeoTransform()
ct = osr.CoordinateTransformation(gcs, pcs)
coordinate = ct.TransformPoint(lon_lat[0], lon_lat[1])
# 地理坐标系转投影坐标系
p_x = np.array([[extend[1], extend[2]], [extend[4], extend[5]]])
p_y = np.array([coordinate[0] - extend[0], coordinate[1] - extend[3]])
# 投影坐标系转图像坐标系
row_col = np.linalg.solve(p_x, p_y)
row = int(float(row_col[1]))
col = int(float(row_col[0]))
img_value = img[row, col]
return img_value
时间匹配
def get_data(_modis_path, xls_path, loc):
"""
:param _modis_path: 图像文件夹地址
:param xls_path: LAI文件地址
:param loc: 要获取的经纬度
:return:
"""
# 获取img文件列表
img_path = []
for i in os.listdir(_modis_path):
if path.splitext(i)[1] == '.img':
img_path.append(path.join(_modis_path, i))
# 获取LAI值和对应时间
data_dict = {'LAI_Value': [],
'LAI_Date': [],
'WDRVI_Path': [],
'WDRVI_Value': []}
wb = xlrd.open_workbook(xls_path)
ws = wb.sheet_by_index(0)
for i in range(0, ws.nrows, 2):
tmp_list1 = ws.row_values(rowx=i, start_colx=3)
tmp_list2 = ws.row_values(rowx=i + 1, start_colx=3)
while '' in tmp_list1:
tmp_list1.remove('')
while '' in tmp_list2:
tmp_list2.remove('')
data_dict['LAI_Value'] += tmp_list1
data_dict['LAI_Date'] += tmp_list2
for m in tmp_list2:
for n in img_path:
x, y = int(m[0:3]), int(n[-17:-14])
if 0 <= (x - y) < 8:
data_dict['WDRVI_Path'].append(n)
data_dict['WDRVI_Value'].append(get_img_value(n, loc))
break
return data_dict
SG滤波及皮尔逊相关系数计算与图表的绘制
if __name__ == "__main__":
modis_path = r"D:\Download\实验一\数据\MODIS"
xls_path = r"D:\Download\实验一\数据\LAI.xls"
loc = [-96.476639, 41.165056]
data = get_data(modis_path, xls_path, loc)
# 利用SG滤波库savgol_filter
x = range(len(data['LAI_Date']))
y1 = data['WDRVI_Value']
y2 = savgol_filter(y1, window_length=7, polyorder=2)
plt.plot(x, y1, color='y', label='pre_filter data')
plt.plot(x, y2, "b--", label="sg result")
plt.legend(loc='lower right')
plt.title('Result of SG Filter')
plt.show()
# 皮尔逊相关系数检验
r, p_value = stats.pearsonr(y1, y2)
print('相关系数为{:.3f},p值为{:.5f}'.format(r, p_value))
# 绘制基本统计量的图表
sg_WDRVI = np.array(y2).reshape(-1, 1)
LAI_value = np.array(data['LAI_Value'])
LAI_value_r = LAI_value.reshape(-1, 1)
data_array = np.hstack([sg_WDRVI, LAI_value_r])
data_frame = pd.DataFrame(data_array, index=None, columns=['WDRVI', 'LAI'])
sns.set(style="ticks", color_codes=True)
plt.figure(dpi=500)
# 绘制WDRVI和LAI的散点图与直方图
sns.pairplot(data=data_frame, vars=['WDRVI', 'LAI'])
plt.show()
plt.savefig('fig.png') # 绘图结果存到本地
# 皮尔逊相关系数矩阵
cmap = sns.diverging_palette(230, 20, as_cmap=True)
sns.heatmap(data_frame.corr(method='pearson'), cmap=cmap, vmax=1, center=0.95,
square=True, linewidths=.5, cbar_kws={"shrink": .5})
plt.title('Pearson Heatmap')
plt.show()
plt.savefig('fig1.png') # 绘图结果存到本地
相关性分析与回归分析
# 一元线性回归
model = LinearRegression()
model.fit(sg_WDRVI, LAI_value)
print('coefficient of determination: {}\n'
'intercept: {}\n'
'slope: {}'.format(model.score(sg_WDRVI, LAI_value), model.intercept_, model.coef_[0]))
y = model.coef_[0] * sg_WDRVI + model.intercept_
plt.scatter(sg_WDRVI, LAI_value, label='raw data', color='r')
plt.legend(loc='lower right')
plt.plot(sg_WDRVI, y, label='fitting', color='b')
plt.legend(loc='lower right')
plt.title('Result of Linear Regression')
plt.show()