Python遥感开发之数据趋势分析Sen+mk

等待着冬天的风

已于 2024-09-07 09:37:54 修改

阅读量2.8k

点赞数 22

CC 4.0 BY-SA版权

分类专栏： Python图像遥感文章标签： Python遥感开发 Python实现Sen+mk

于 2024-05-19 22:23:46 首次发布

本文链接：https://blog.csdn.net/qq_32306361/article/details/139049555

Python遥感开发之数据趋势分析Sen+mk

1 方法介绍
- 1.1 Theil-Sen Median方法
- 1.2 Mann-Kendall方法
2 Python实现Sen+mk
- 2.1 代码（作者实现，推荐使用）
- 2.2 另外版本
3 最终结果

前言：主要使用Python完成遥感时间序列数据趋势分析Sen+mk，得到slope、trend、p、s、tau、z指标。

1 方法介绍

1.1 Theil-Sen Median方法

又被称为 Sen 斜率估计，是一种稳健的非参数统计的趋势计算方法。该方法计算效率高，对于测量误差和离群数据不敏感，常被用于长时间序列数据的趋势分析中。对于后续代码计算结果中的slope.tif解读,当slope大于0表示随时间序列呈现上升趋势；slope小于0表示随时间序列呈现下降趋势。
在这里插入图片描述

1.2 Mann-Kendall方法

Mann-Kendall是一种非参数统计检验方法，最初由Mann在1945年提出，后由Kendall和Sneyers进一步完善，其优点是不需要测量值服从正态分布，也不要求趋势是线性的，并且不受缺失值和异常值的影响，在长时间序列数据的趋势显著检验中得到了十分广泛的应用。对于后续代码计算结果中的z.tif,当|Z|大于1.65、1.96和2.58时，表示趋势分别通过了置信度为90%、95%和99%的显著性检验。
在这里插入图片描述

2 Python实现Sen+mk

2.1 代码（作者实现，推荐使用）

# coding:utf-8
import numpy as np
import pymannkendall as mk
import os
from osgeo import gdal


def read_tif(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize  # 图像长度
    row = dataset.RasterYSize  # 图像宽度
    geotrans = dataset.GetGeoTransform()  # 读取仿射变换
    proj = dataset.GetProjection()  # 读取投影
    data = dataset.ReadAsArray()  # 转为numpy格式
    data = data.astype(np.float32)
    no_data_value = data[0][0]  # 设定NoData值
    data[data == no_data_value] = np.nan  # 将NoData转换为NaN
    return col, row, geotrans, proj, data


def save_tif(data, reference_file, output_file):
    ds = gdal.Open(reference_file)
    shape = data.shape
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(output_file, shape[1], shape[0], 1, gdal.GDT_Float32)  # 保存的数据类型
    dataset.SetGeoTransform(ds.GetGeoTransform())
    dataset.SetProjection(ds.GetProjection())
    dataset.GetRasterBand(1).WriteArray(data)
    dataset.FlushCache()


def get_tif_list(directory):
    """
    获取目录下所有TIF文件的完整路径
    """
    return [os.path.join(directory, f) for f in os.listdir(directory) if f.endswith('.tif')]


if __name__ == '__main__':
    # 输入文件夹路径
    input_directory = r"E:\AAWORK\学校相关\rsei\data\四季\rsei-winter"
    output_directory = r"E:\AAWORK\学校相关\rsei\data\四季\sen+mk"

    # 读取文件列表
    file_list = get_tif_list(input_directory)

    # 读取第一个文件获取基本信息（列、行、仿射变换、投影）
    col, row, geotrans, proj, _ = read_tif(file_list[0])

    # 预分配结果数组
    slope_data = np.zeros((row, col), dtype=np.float32)
    s_data = np.zeros((row, col), dtype=np.float32)
    z_data = np.zeros((row, col