【基于Python的栅格数据逐像元Sen+mk趋势分析】实现一键运行出结果，简单易行

前言

结合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值。

注意事项：

1. 数据准备： 确保输入的遥感数据格式正确，输入数据应是使用Arcmap波段合成工具后包含多个时间点的多波段栅格数据。如多年 NDVI 数据，每一年是一个波段。

2. 代码修改： 在本代码中仅修改输入输出路径即可运行出结果（第56-61行，代码中有明确注释）。

3. 结果解释： 了解趋势分析结果的含义，特别是Sen's slope和显著性z-value。这有助于准确解释分析的实际意义。

4. 空间分辨率： 在处理栅格数据时，注意空间分辨率的影响。高分辨率数据可能需要更多计算资源，而低分辨率数据可能无法捕捉到细微的趋势。

5. 库版本： 确保使用的pymannkendall和rasterio库是最新版本，以获得可能的性能和功能改进。

欢迎交流，联系邮箱：83827730@qq.com 。