【Python】遥感数据趋势分析Sen+mk

方法介绍

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

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

代码介绍

此代码大部分来自Github大佬分享，我修改了报错的个别代码后亲自测试pycharm下只需改路径就可以运行：

# coding:utf-8 ''' 已全部实现 ''' import numpy as np import pymannkendall as mk import os import rasterio as ras def sen_mk_test(image_path, outputPath): # image_path:影像的存储路径 # outputPath:结果输出路径 global path1 filepaths = [] for file in os.listdir(path1): filepath1 = os.path.join(path1, file) filepaths.append(filepath1) # 获取影像数量 num_images = len(filepaths) # 读取影像数据 img1 = ras.open(filepaths[0]) # 获取影像的投影，高度和宽度 transform1 = img1.transform height1 = img1.height width1 = img1.width array1 = img1.read() img1.close() # 读取所有影像 for path1 in filepaths[1:]: if path1[-3:] == 'tif': print(path1) img2 = ras.open(path1) array2 = img2.read() array1 = np.vstack((array1, array2)) img2.close() nums, width, height = array1.shape # 写影像 def writeImage(image_save_path, height1, width1, para_array, bandDes, transform1): with ras.open( image_save_path, 'w', driver='GTiff', height=height1, width=width1, count=1, dtype=para_array.dtype, crs='+proj=latlong', transform=transform1, ) as dst: dst.write_band(1, para_array) dst.set_band_description(1, bandDes) del dst # 输出矩阵，无值区用-9999填充 slope_array = np.full([width, height], -9999.0000) z_array = np.full([width, height], -9999.0000) Trend_array = np.full([width, height], -9999.0000) Tau_array = np.full([width, height], -9999.0000) s_array = np.full([width, height], -9999.0000) p_array = np.full([width, height], -9999.0000) # 只有有值的区域才进行mk检验 c1 = np.isnan(array1) sum_array1 = np.sum(c1, axis=0) nan_positions = np.where(sum_array1 == num_images) positions = np.where(sum_array1 != num_images) # 输出总像元数量 print("all the pixel counts are {0}".format(len(positions[0]))) # mk test for i in range(len(positions[0])): print(i) x = positions[0][i] y = positions[1][i] mk_list1 = array1[:, x, y] trend, h, p, z, Tau, s, var_s, slope, intercept = mk.original_test(mk_list1) ''' trend: tells the trend (increasing, decreasing or no trend) h: True (if trend is present) or False (if trend is absence) p: p-value of the significance test z: normalized test statistics Tau: Kendall Tau s: Mann-Kendal's score var_s: Variance S slope: Theil-Sen estimator/slope intercept: intercept of Kendall-Theil Robust Line ''' if trend == "decreasing": trend_value = -1 elif trend == "increasing": trend_value = 1 else: trend_value = 0 slope_array[x, y] = slope # senslope s_array[x, y] = s z_array[x, y] = z Trend_array[x, y] = trend_value p_array[x, y] = p Tau_array[x, y] = Tau all_array = [slope_array, Trend_array, p_array, s_array, Tau_array, z_array] slope_save_path = os.path.join(result_path, "slope.tif") Trend_save_path = os.path.join(result_path, "Trend.tif") p_save_path = os.path.join(result_path, "p.tif") s_save_path = os.path.join(result_path, "s.tif") tau_save_path = os.path.join(result_path, "tau.tif") z_save_path = os.path.join(result_path, "z.tif") image_save_paths = [slope_save_path, Trend_save_path, p_save_path, s_save_path, tau_save_path, z_save_path] band_Des = ['slope', 'trend', 'p_value', 'score', 'tau', 'z_value'] for i in range(len(all_array)): writeImage(image_save_paths[i], height1, width1, all_array[i], band_Des[i], transform1) # 调用 path1 = r"E:\Test\ACMEI\ACMEI" result_path = r"E:\Test\ACMEI\ACMEImk" sen_mk_test(path1, result_path)

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• 52
• 53
• 54
• 55
• 56
• 57
• 58
• 59
• 60
• 61
• 62
• 63
• 64
• 65
• 66
• 67
• 68
• 69
• 70
• 71
• 72
• 73
• 74
• 75
• 76
• 77
• 78
• 79
• 80
• 81
• 82
• 83
• 84
• 85
• 86
• 87
• 88
• 89
• 90
• 91
• 92
• 93
• 94
• 95
• 96
• 97
• 98
• 99
• 100
• 101
• 102
• 103
• 104
• 105
• 106
• 107
• 108
• 109
• 110
• 111
• 112
• 113
• 114
• 115
• 116
• 117
• 118
• 119
• 120
• 121
• 122

后续处理三步走

1.对于生成的tif栅格，只需要取slope.tif以及z.tif进行重分类。

slope>0赋值1表示增加 |z|>1.96赋值2表示显著

slope=0赋值0表示不变 |z|<=1.96赋值1表示不显著

slope<0赋值-1表示减少

（此处还可根据显著性阈值再细分为更多类，我仅作基本演示分为5类）

2.再相乘得到

-2:显著减少;-1:不显著减少;0:稳定不变;1:不显著增加; 2:显著增加

3.在Arcmap出图即可！