基于Python(Arcpy)对遥感影像（TIF）的Mann-Kendall算法突变检测实现

主体函数：

由多张影像先求每年的平均值，每年对应一个平均值，对这些值进行突变检测。

#coding=utf-8
import time
start = time.clock()
import arcpy
import math
import matplotlib.pyplot as plt
# 对每年的栅格数据求一个平均数值
def mean_year(workspace,zonal_mask):
    arcpy.env.workspace = workspace
    arcpy.env.overwriteOutput = True
    arcpy.CheckOutExtension("spatial")
    raster_list = arcpy.ListRasters("*","tif")
    name_table = "Temp_table"
    mean_ls = []
    if arcpy.Exists(name_table):
        arcpy.Delete_management(name_table)
    for raster in raster_list:
        #print raster
        arcpy.sa.ZonalStatisticsAsTable(in_zone_data=zonal_mask,\
                                        zone_field="name",\
                                        in_value_raster=raster,\
                                        out_table=name_table,\
                                        statistics_type="MEAN")
        cursor = arcpy.da.SearchCursor(name_table,["MEAN"])
        for row1 in cursor:
            mean = row1[0]
            #print mean
            mean_ls.append(mean)
    return mean_ls
# 计算UFK与UBK序列
def cal_UFKUBK(mean_ls):
    # 计算UFK 正序列计算完成
    s = 0
    # Sk = []
    UFk = []
    UFk.append(0)
    n = len(mean_ls)
    for i in range(1,n,1):
        for j in range(0,i,1):
            if mean_ls[i] > mean_ls[j]:
                s = s +1
            else:
                s = s + 0
        # Sk.append(s)
        ESk = (i+1)*(i)/4.0
        Var = (i+1)*i*(2*(i+1)+5)/72.0
        UF = (s-ESk)*1.0/math.sqrt(Var)
        UFk.append(UF)
    print UFk
    # 计算UBK 逆序列计算完成
    mean_lsT = list(reversed(mean_ls))
    s_T = 0
    UBk_T = []
    UBk_T.append(0)
    for i in range(1,n,1):
        for j in range(0,i,1):
            if mean_lsT[i] > mean_lsT[j]:
                s_T = s_T + 1
            else:
                s_T = s_T + 0
        ESk_T = (i + 1) * (i) / 4.0
        Var_T = (i + 1) * i * (2 * (i + 1) + 5) / 72.0
        UB_T = (s_T - ESk_T) * 1.0 / math.sqrt(Var_T)
        UBk_T.append(-1*UB_T)

    UBk = list(reversed(UBk_T))
    res_ls = []
    res_ls.append(UFk)
    res_ls.append(UBk)
    return res_ls
def cal_changePoint(res_ls,year_ls):
    # 找出突变点
    D = []
    change_point = []
    n = len(year_ls)
    for i in range(n):
        D.append(res_ls[0][i] - res_ls[1][i])
    for i in range(1, n, 1):
        if D[i - 1] * D[i] < 0:
            if math.fabs(D[i - 1]) > math.fabs(D[i]):
                change_point.append(i+year_ls[0])
            else:
                change_point.append(i - 1+year_ls[0])
    point = list(set(change_point))
    print "突变点为："
    print point
    # 画出UFK与UBK的趋势图
    plt.plot(year_ls, res_ls[0], label="UFk", color="red")
    plt.plot(year_ls, res_ls[1], label="UBk", color="blue")

    x_lim = plt.xlim()
    plt.plot(x_lim, [-1.96, -1.96], "m--", color='r')
    plt.plot(x_lim, [0, 0], "m--")
    plt.plot(x_lim, [1.96, 1.96], "m--", color='r')

    plt.legend(loc="upper right")
    plt.ylabel("UFk-UBk")
    plt.show()
# 主函数
if __name__ == "__main__":
    path_tif = "F:/GPP_Temp/res_mean"
    zonal_mask = "F:/shp_water/china_water.shp"
    year_ls = []
    for year in range(1982,2016,1):
        year_ls.append(year)
    dt = mean_year(path_tif,zonal_mask)
    res_ls = cal_UFKUBK(dt)
    cal_changePoint(res_ls,year_ls)

总结：使用时把主函数中的数路径与年序列改成自己对应的。