[人工智能]python实现MK趋势分析（不借助MK的第三方库）

python中有MK趋势分析的第三方库，但我觉得有时候也应该自己实现这种算法，个人觉得这个代码符合算法逻辑，如果各位看官觉得有问题请与我联系！谢谢！

# -*- coding: utf-8 -*-
"""
@Time    : 2022/3/22 13:01
@Author  : FJC
@File    : MK趋势分析.py
@Software: win10  python3.7
"""
import numpy as np
import os
from osgeo import gdal

path = r"G:\result\threshold\SOS_TIF"


def m_k(images, out_image):
    """

    :param images: 输入影像
    :param out_image:输出影像
    """
    open_tif_1 = gdal.Open(r"G:\result\threshold\SOS_TIF\2000.tif")
    rows = open_tif_1.RasterYSize  # 行数
    cols = open_tif_1.RasterXSize  # 列数
    images_pixels = np.zeros((21, rows, cols))
    k = 0
    for image in images:
        image = str(image)
        open_tif_2 = gdal.Open(image)
        band = open_tif_2.GetRasterBand(1).ReadAsArray()  # 获取波段的矩阵
        images_pixels[k, :, :] = band
        k += 1
    data = np.zeros((210, rows, cols))
    n = 0
    for i in range(20):
        for j in range(i+1,21):
            data[n, :, :] = images_pixels[j, :, :] - images_pixels[i, :, :]
            n += 1
    sgn = np.where(data < 0, -1, np.where(data > 0, 1, 0))
    s = np.sum(sgn, axis=0)
    var = 33.116
    zc = np.where(s > 0, (s - 1) / var, np.where(s < 0, (s + 1) / var, 0))
    driver = gdal.GetDriverByName("GTiff")
    creat_tif = driver.Create(out_image, cols, rows, 1, gdal.GDT_Float32)
    creat_tif.SetProjection(open_tif_1.GetProjection())
    creat_tif.SetGeoTransform(open_tif_1.GetGeoTransform())
    out_tif = creat_tif.GetRasterBand(1)
    write_band = out_tif.WriteArray(zc)
    print("M_K趋势分析计算完成")


if __name__ == "__main__":
    out_image = r"G:\result\threshold\SOS_result\SOS_MK.tif"  # 定义一个栅格名
    image_names = os.listdir(path)
    images_paths = []
    for images_name in image_names:
        images_path = path + '//' + images_name
        images_paths.append(images_path)
    print(images_paths)
    m_k(images_paths, out_image)

最后我发现这个数据量太大了，应该还需改进，数据量小的话应该是可以的。