Python实现时间序列栅格数据MK检验(Mann-Kendall)

01.Mann-Kendall 趋势检验

Mann-Kendall 趋势检验（有时称为 MK 检验）用于分析时间序列数据，以确定持续增加或减少的趋势（单调趋势）。它是一种非参数检验，这意味着它适用于所有分布（即数据不必满足正态性假设），但数据不应具有序列相关性。如果数据具有序列相关性，则可能会影响显著水平（p 值）。

02.pyMannKendall

pyMannKendall是用于非参数 Mann-Kendall 趋势测试系列的 python 包，本文中使用它实现的MK检验。一维数组的MK检验代码如下所示：

import pymannkendall as mk
import numpy as np

#创建数据
data = np.array([34, 33, 29,31,30,39,27, 26, 27, 25, 23, 21, 17, 19, 16, 14, 16, 12, 11, 9, 8])
# MK检验
mk.original_test(data)

# 输出结果
# Mann_Kendall_Test(
#    trend='decreasing', 
#    h=True, 
#    p=1.5860952995438993e-08, 
#    z=-5.651980563258717, 
#    Tau=-0.8952380952380953, 
#    s=-188.0, 
#    var_s=1094.6666666666667, 
#    slope=-1.3333333333333333, 
#    intercept=36.33333333333333)

结果解释：

trend：趋势增加、减少或无趋势（increasing, decreasing or no trend）

h：True（趋势存在）或 False（趋势不存在）

p：显著性检验的p值

z：归一化检验统计量

Tau：Kendall Tau

s：Mann-Kendall分数

var_s：方差 S

slop：Theil-Sen 斜率

intercept：截距

03.时间序列栅格数据MK检验

import glob
import os
import numpy as np
import rasterio as rio
import pymannkendall as mk

# MK检验
def apply_mk_func(arr):
    if len(arr) > 1 and not np.any(np.isnan(arr)):
        trend, h, p, z, Tau, s, var_s, slope, intercept = mk.original_test(arr)
        # 将趋势increasing, decreasing , no trend，替换为 1、-1、0
        trend_value = 0
        if trend == "decreasing":
            trend_value = -1
        elif trend == "increasing":
            trend_value = 1
        return [trend_value, h, p, z, Tau, s, var_s, slope, intercept]
    return np.full((9),np.nan)

# 输出栅格数据
def write_image(image_save_path, data_array, out_meta):
    with rio.open(
            image_save_path,
            'w',
           **out_meta,
    ) as dst:
        dst.write_band(1, data_array)
    del dst

if __name__ == "__main__":
    
    in_path = "data\\ndvi"
    out_path = "data\\ndvi_mk"
    
    file_list =  glob.glob(in_path+'\\*.tif')
    count = len(file_list)
    
    ds_f = rio.open(file_list[0])
    # 获取mate,后续导出tif时候使用；
    meta = ds_f.meta
    
    width = ds_f.width
    height = ds_f.height
    # 构建空三维数组，用于组建时间序列栅格数据；
    time_array = np.zeros((count, height, width))
    for i,file in enumerate(file_list):
        print("读取数据：{0}：{1}".format(i+1, file))
        ds = rio.open(file)
        data = ds.read(1)
        mask = ds.read_masks(1)
        data[mask == 0] = np.nan
        time_array[i] = data
    # 时间轴向应用mk检验函数
    print("=====开始mk检验=====")
    mk_result = np.apply_along_axis(apply_mk_func, 0, time_array)
    out_list = ["trend_value", "h", "p", "z", "Tau", "s", "var_s", "slope", "intercept"]
    print('=====开始导出数据=====')
    for i, out_file in enumerate(out_list):
        image_save_path = os.path.join(out_path, "{}.tif".format(out_file) )
        data_array = mk_result[i]
        write_image(image_save_path, data_array, meta)
        print("导出'{0}'数据：{1}".format(out_file, image_save_path))
    print("完成！！！")

参考文章

[1].趋势检验方法（二）MK趋势检验_mk检验-CSDN博客

[2].MK+Sen趋势检验（长时间栅格数据）_theil-sen median图像行列不一致-CSDN博客

[3].pymannkendall · PyPI