前言
结合Sens斜率和Mann-Kendall趋势检测(MK)方法,使用pymannkendall
和rasterio
库。这两者的结合为我们提供了一种强大的工具,能够深入挖掘遥感数据中的趋势和显著性。
pymannkendall
提供了Sen+mk趋势检测的Python实现。这个库背后的算法经过数学验证,是Mann-Kendall检测的标准方法之一。在科学研究和数据分析社区中,pymannkendall
被广泛采用,许多研究项目都选择使用它来进行趋势分析。
rasterio
是一个用于处理栅格数据的强大库,广泛应用于遥感和地理信息系统领域。该库提供了高效的栅格数据读取和写入功能,同时支持处理空间信息和坐标参考系统。rasterio
的可靠性得益于其在开源社区中的长期维护和不断改进。
话不多说,先上代码:
代码:
import pymannkendall as mk
import numpy as np
import rasterio
def read_raster(file_path):
"""
读取栅格数据
Parameters:
- file_path: 栅格数据文件路径
Returns:
- data: 读取到的栅格数据
"""
with rasterio.open(file_path) as src:
return src.read()
def sen_mk_map(data):
"""
计算Sen'slope和显著性 z_value 的栅格数据
Parameters:
- data: 包含多个波段的栅格数据,例如 2000-2020多年NDVI,使用Arcmap中波段合成工具,将21年NDVI数据合称为具有21个波段的单栅格数据。
Returns:
- slope_value_map: 存储slope值的栅格数据
- z_value_map: 存储显著性 z-value 的栅格数据
"""
num_years, rows, cols = data.shape
# 初始化存储结果的数组
slope_value_map = np.zeros((rows, cols), dtype=np.float32)
z_value_map = np.zeros((rows, cols), dtype=np.float32)
# 遍历每个像元位置
for i in range(rows):
for j in range(cols):
# 提取对应位置的值
values_data = data[:, i, j]
# 计算slope和显著性 z-value
trend, h, p, z, Tau, s, var_s, slope, intercept = mk.original_test(values_data)
# 获取slope和显著性 z_value
slope_value = slope
z_value = z
# 将计算结果存储在对应位置
slope_value_map[i, j] = slope_value
z_value_map[i, j] = z_value
return slope_value_map, z_value_map
# 输入数据路径
ndvi_path = "E:\path"
# slope栅格输出数据路径
slope_path="E:\path"
# z栅格输出数据路径
z_path="E:\path"
# 读取栅格数据
data_stack = read_raster(ndvi_path)# 多年 NDVI 数据,每一年是一个波段
with rasterio.open(ndvi_path) as src:
# 计算slope和显著性 z-value 栅格数据
slope_value_map, z_value_map = sen_mk_map(data_stack)
# 将结果写入新的栅格数据文件 - 存储slope值
output_path_corr = slope_path
with rasterio.open(
output_path_corr,
'w',
driver='GTiff',
height=src.height,
width=src.width,
count=1,
dtype=src.dtypes[0], # 使用输入文件的数据类型
crs=src.crs,
transform=src.transform,
) as dst_slope:
dst_slope.write(slope_value_map, 1)
# 将显著性 z-value 结果写入新的栅格数据文件
output_path_z_value = z_path
with rasterio.open(
output_path_z_value,
'w',
driver='GTiff',
height=src.height,
width=src.width,
count=1,
dtype=src.dtypes[0], # 使用输入文件的数据类型
crs=src.crs,
transform=src.transform,
) as dst_z_value:
dst_z_value.write(z_value_map, 1)
遥感数据读取和处理:
使用rasterio
库读取遥感数据,该库提供了强大的功能,包括读取、写入、处理空间信息等。在这里,通过自定义read_raster
函数读取了栅格数据,然后通过自定义函数sen_mk_map计算Sen'slope和显著性z-value的栅格数据。
Sens斜率计算:
Sens斜率是一种用于衡量时间序列数据中趋势的方法。通过计算Sens斜率,我们可以了解到遥感数据在时间上的变化趋势。使用pymannkendall
库的mk.original_test函数来计算Sens斜率。
Mann-Kendall趋势检测:
Mann-Kendall趋势检测是一种非参数统计方法,用于检测时间序列中的趋势。通过比较数据点的大小关系,该方法可以判断数据是否呈现出显著的趋势。使用pymannkendall
库mk.original_test函数返回z值。
注意事项:
-
数据准备: 确保输入的遥感数据格式正确,输入数据应是使用Arcmap波段合成工具后包含多个时间点的多波段栅格数据。如多年 NDVI 数据,每一年是一个波段。
-
代码修改: 在本代码中仅修改输入输出路径即可运行出结果(第56-61行,代码中有明确注释)。
-
结果解释: 了解趋势分析结果的含义,特别是Sen's slope和显著性z-value。这有助于准确解释分析的实际意义。
-
空间分辨率: 在处理栅格数据时,注意空间分辨率的影响。高分辨率数据可能需要更多计算资源,而低分辨率数据可能无法捕捉到细微的趋势。
-
库版本: 确保使用的
pymannkendall
和rasterio
库是最新版本,以获得可能的性能和功能改进。
欢迎交流,联系邮箱:83827730@qq.com 。