Python | 栅格数据Sen+MK长时间序列趋势分析+显著性检验保姆级教程_theil-sen median趋势分析

 遥感代码星球的第002篇代码分享 

今天给大家分享的是:

Theil-Sen Median趋势分析

+

Mann-Kendall显著性检验

基于遥感数据的植被碳水利用效率时空变化和归因分析[J].生态学报,2024,44(01):377-391.DOI:10.20103/j.stxb.202207232112.

在研究气候变化、环境变化、生态变化等地方时，如何准确分析时间序列数据的趋势是一个重要的课题。Theil-Sen Median趋势分析（简称Sen分析）和Mann-Kendall显著性检验（简称MK检验）是两种常用的非参数方法，它们在处理气候、环境、生态等地方的趋势分析中有着广泛的应用。本文将详细介绍这两种方法的基本原理，并探讨它们如何结合使用以提高趋势分析的准确性。同时结尾附完整代码及案例数据供大家学习使用。

01 原理及公式

1. Theil-Sen Median 趋势分析（Sen分析）

Theil-Sen Median 趋势分析是一种稳健的非参数统计方法，适用于具有异常值或缺失值的时间序列数据。该方法通过计算数据的中位数斜率（即Sen斜率）来评估时间序列的趋势。不同于传统的最小二乘法，Sen分析通过中位数的方式计算斜率，使其对异常值和极端值不敏感，从而能够更稳健地捕捉数据中的长期趋势。

具体的计算方法为(本文以ET为例)：

其中，ETi和ETj分别为时间序列中不同时间点的数据值Sen为中位数斜率。

若Sen>0，则表示数据存在上升趋势；若Sen<0，则表示数据存在下降趋势。

2. Mann-Kendall 显著性检验（MK检验）

Mann-Kendall检验是一种用于检测时间序列趋势的非参数统计方法。与Sen分析不同，MK检验的目的是判断时间序列的趋势是否显著，它不要求数据服从任何特定的分布。MK检验通过对时间序列中每对数据点进行比较，计算出趋势的符号，并利用统计量评估数据的趋势方向和显著性。

MK检验的基本步骤如下：

1. 对时间序列中所有的每对数据点进行比较，得到每对数据的符号（升序为1，降序为-1，相等为0）；

2. 计算符号函数的累加和，得到统计量S：

    

   

   其中，sgn(ETj−ETi)为符号函数。

3. MK检验的检验统计量S反映了时间序列的总体趋势。根据S值的大小及其方差，可以进一步计算出Z值，用于判断ET变化的显著性。计算S的方差，进一步得到Z值，用于判断趋势的显著性：

    

   

   Z值越大，表示趋势越显著。

    

3. MK显著性检验结果的划分

在进行MK检验时，Z值的大小决定了趋势的显著性水平。在常见的0.05置信水平下，MK检验的显著性结果可以通过Z值来判断：

• 显著变化: 上升趋势（Z>1.96）或下降趋势（Z<−1.96）

• 不显著变化：（−1.96≤Z≤1.96）

4. Sen+MK方法的优势与应用

• 稳健性强：Sen分析能够有效处理含有异常值或缺失值的时间序列数据，适用于多种复杂环境下的数据分析。

• 显著性检验：MK检验能够评估趋势的显著性，确保得到的趋势结果是有统计学意义的。

• 适用范围广泛：Sen+MK方法在气候变化、环境监测、生态变化等地方中广泛应用，尤其适用于无法满足正态分布假设的时间序列数据。

02 代码复现

1. 代码可实现：

✅逐像元栅格进度可视化

✅只需要改动标注有✏️的位置

✅可1比1复刻

2. 文件夹结构：

3.具体代码：

import osimport rasterioimport numpy as npfrom tqdm import tqdmimport matplotlib.pyplot as pltfrom pathlib import Path# 获取当前工作目录base_path = os.getcwd()# 定义数据和结果路径data_path = os.path.join(base_path, \'data\') # ✏️result_path = os.path.join(base_path, \'results\') # ✏️# 创建结果目录（如果不存在）os.makedirs(result_path, exist_ok=True)# 定义年份范围start_year = 2000 # ✏️end_year = 2020 # ✏️cd = end_year - start_year + 1 # 时间跨度try: # 读取第一个栅格文件以获取元数据和尺寸信息 first_file = os.path.join(data_path, f\'{start_year}.tif\') with rasterio.open(first_file) as src: a = src.read(1) # 读取栅格数据 transform = src.transform # 栅格的空间转换信息 metadata = src.meta.copy() # 获取元数据，包括GeoTIFFTags # 修复：更新元数据中的数据类型为float32，解决dtype不匹配问题 metadata.update({ \'dtype\': \'float32\', \'nodata\': np.nan # 设置空值为NaN }) # 获取栅格尺寸 m, n = a.shape print(f\"栅格尺寸: {m} x {n}\") # 创建空数组来存储所有年份的数据 # 使用三维数组以避免展平和重塑操作 all_data = np.full((m, n, cd), np.nan) # 加载每一年的数据 print(\"正在加载每一年的栅格数据...\") for i, year in enumerate(tqdm(range(start_year, end_year + 1), desc=\"加载年度数据\")): filename = os.path.join(data_path, f\'{year}.tif\') with rasterio.open(filename) as src: all_data[:, :, i] = src.read(1) # 创建结果数组 sen_result = np.full((m, n), np.nan) # 使用Sen\'s Slope方法计算趋势 # 这是计算最耗时的部分，使用tqdm显示进度 print(\"正在计算Sen\'s Slope趋势...\") valid_pixels = 0 # 优化：仅处理有效像素（至少有一个非空值的像素） valid_mask = np.any(all_data > 0, axis=2) total_valid = np.sum(valid_mask) for i in tqdm(range(m), desc=\"计算Sen\'s Slope\"): for j in range(n): if valid_mask[i, j]: data = all_data[i, j, :] if np.all(~np.isnan(data)) and np.min(data) > 0: # 检查数据是否包含有效值  slopes = []  for k1 in range(1, cd): for k2 in range(k1): # 计算变化率 slope = (data[k1] - data[k2]) / (k1 - k2) slopes.append(slope)  # 取中位数作为最终斜率值  if slopes: # 确保slopes不为空 sen_result[i, j] = np.median(slopes) valid_pixels += 1 print(f\"有效像素数: {valid_pixels}/{total_valid} (共{m*n}像素)\") # 设置输出路径并将Sen\'s Slope结果保存为GeoTIFF sen_output_path = os.path.join(result_path, \'基于sen的ET变化趋势.tif\') # ✏️ # 确保结果为float32类型 sen_result = sen_result.astype(\'float32\') # 保存Sen\'s Slope结果 with rasterio.open(sen_output_path, \'w\', **metadata) as dst: dst.write(sen_result, 1) print(\'Sen\\\'s Slope处理完成! 结果已保存至:\', sen_output_path) # 计算Mann-Kendall检验 print(\"正在计算Mann-Kendall检验...\") mk_result = np.full((m, n), np.nan) for i in tqdm(range(m), desc=\"计算Mann-Kendall检验\"): for j in range(n): if valid_mask[i, j]: data = all_data[i, j, :] if np.all(~np.isnan(data)) and np.min(data) > 0: # 检查数据是否包含有效值  sgnsum = []  for k1 in range(1, cd): for k2 in range(k1): # 计算符号差异 sgn = np.sign(data[k1] - data[k2]) sgnsum.append(sgn)  # 计算符号差异的总和  mk_result[i, j] = np.sum(sgnsum) # 计算Z值 print(\"正在计算Z值...\") vars_mk = cd * (cd - 1) * (2 * cd + 5) / 18 z_scores = np.full((m, n), np.nan) # 处理不同情况的Z值计算 z_scores[~np.isnan(mk_result) & (mk_result == 0)] = 0 z_scores[~np.isnan(mk_result) & (mk_result > 0)] = (mk_result[~np.isnan(mk_result) & (mk_result > 0)] - 1) / np.sqrt(vars_mk) z_scores[~np.isnan(mk_result) & (mk_result < 0)] = (mk_result[~np.isnan(mk_result) & (mk_result < 0)] + 1) / np.sqrt(vars_mk) # 保存Mann-Kendall检验结果 mk_output_path = os.path.join(result_path, \'MK检验结果.tif\') # ✏️ z_scores = z_scores.astype(\'float32\') with rasterio.open(mk_output_path, \'w\', **metadata) as dst: dst.write(z_scores, 1) print(\'Mann-Kendall检验处理完成! 结果已保存至:\', mk_output_path) # 对Sen\'s slope和MK检验结果进行重分类 print(\"正在进行重分类...\") # 创建重分类结果数组 S2 = np.full((m, n), np.nan) M2 = np.full((m, n), np.nan) # 重分类Sen\'s slope结果 S2[np.isnan(sen_result)] = -9999 S2[~np.isnan(sen_result) & (sen_result = 0.0005)] = 1 S2[~np.isnan(sen_result) & (sen_result > -0.0005) & (sen_result  1.96] = 2 M2[~np.isnan(z_scores) & (z_scores <= 1.96)] = 1 # 计算最终重分类结果 reclassify = (S2 * M2).astype(np.int16) # 设置输出路径并将重分类结果保存为GeoTIFF reclass_output_path = os.path.join(result_path, \'重分类结果.tif\') # ✏️ # 更新元数据中的数据类型为int16 metadata.update({ \'dtype\': \'int16\', \'nodata\': -9999 }) # 保存重分类结果 with rasterio.open(reclass_output_path, \'w\', **metadata) as dst: dst.write(reclassify, 1) print(\'重分类处理完成! 结果已保存至:\', reclass_output_path) except Exception as e: print(f\"处理过程中发生错误: {str(e)}\") 

4. 可视化进度

5.输出结果

6. 本文示例数据：

MOD16A2H 2000-2020年宁夏部分

本文可视化部分的代码以 Jupyter Notebook的

.ipynb文件展现

你也可以将其中内容复制到新的.py文件中

本人代码能力有限

欢迎同学们指出错误或优化意见

我们一起成长进步

❤

 关注遥感代码星球 

~评论区回复 C002 获取代码及示例数据链接~