基于python在算法流程改进的MK趋势分析（不借助第三方库）

刚刚上一篇由于不能运行大数据，导致用途有限，但算法上是没问题的。这个程序是对上篇算法上的改进。个人觉得算法上符合MK算法，如果各位看官发现bug请指正，谢谢！

# -*- 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:resultthresholdSOS_TIF"


def m_k(images, out_image):
    """

    :param images: 输入影像
    :param out_image:输出影像
    """
    open_tif_1 = gdal.Open(r"G:resultthresholdSOS_TIF2000.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
    s_1 = np.zeros((20, rows, cols))
    n = 20
    for i in range(20):
    	k = 0
        data = np.zeros((n, rows, cols))
        for j in range(i + 1, 21):
            data[k, :, :] = images_pixels[j, :, :] - images_pixels[i, :, :]
            k += 1
        sgn = np.where(data < 0, -1, np.where(data > 0, 1, 0))
        s_1[i, :, :] = np.sum(sgn, axis=0)  # 先进行第一层累加求和
        n -= 1
    s = np.sum(s_1, 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:resultthresholdSOS_resultSOS_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)