基于Python(Arcpy)批量对遥感影像MK与SEN检验，结果为两张TIF图（MK与SEN）

def mkSen_cal(source_folder,outpt_folder):
    arcpy.CheckOutExtension("spatial")
    arcpy.env.workspace = source_folder
    print "源数据的路径为："+str(arcpy.env.workspace)
    arcpy.overwriteOutput = True
    arcpy.CheckOutExtension("spatial")
    muti_path = outpt_folder
    print "成果保存的路径为:"+muti_path
    raster_list = arcpy.ListRasters("*","tif")
    print "tif文件的个数："+str(len(raster_list))
    print raster_list
    sum_array = 0
    n = len(raster_list)

    print "开始求统计量S"
    for j in range(n-1):
        for k in range(j+1,n,1):
        #print str(k)+"--"+str(j)
            raster1 = arcpy.env.workspace + "/"+raster_list[j]
            raster2 = arcpy.env.workspace + "/"+raster_list[k]
            desc1 = arcpy.Describe(raster1)
            desc2 = arcpy.Describe(raster2)
            nameRaster = "正在叠加:"+"{0}_{1}.tif".format(desc1.basename,desc2.basename)
            print nameRaster

            arr = Muti_2page.muti_tif(raster1,raster2)
            sum_array = sum_array + arr
    arcpy.env.workspace = muti_path
    s_Raster_path = "s_Raster.tif"
    print s_Raster_path+"统计量计算成功"
    if arcpy.Exists(s_Raster_path):
        arcpy.Delete_management(s_Raster_path)
        print "正在删除文件......"
    sum_array.save(s_Raster_path)
    print "S统计量栅格计算并保存成功"

    print "开始计算Z统计量"
    s_Raster = arcpy.Raster(muti_path+"/"+"s_Raster.tif")

    var = math.sqrt(n*(n-1)*(2*n+5)/18)
    array1 = arcpy.sa.Con(s_Raster,(s_Raster+1)/var,s_Raster,"VALUE < 0")
    array2 = arcpy.sa.Con(array1,(s_Raster-1)/var,array1,"VALUE > 0")
    z_Raster_path = muti_path+"/"+"z_Raster.tif"
    arcpy.Delete_management(array1)
    if arcpy.Exists(z_Raster_path):
        arcpy.Delete_management(z_Raster_path)
    array2.save(z_Raster_path)
    print "计算Z统计量成功"
    print "开始对Z统计量进行评估"
    Muti_2page.z_belive(z_Raster_path, muti_path)
    print "Z统计量的评估完成"

    print "Sen计算倾斜度"
    slope1 = []
    slope_raster_path = muti_path + "/" + "slope_raster.tif"
    print  "n为："+str(n)
    for k in range(0, n - 1):
        for j in range(k + 1, n):
            raster_j = source_folder + "/" + raster_list[j]
            raster_k = source_folder + "/" +raster_list[k]
            desc1 = arcpy.Describe(raster_j)
            desc2 = arcpy.Describe(raster_k)
            nameRaster = "正在叠加:" + "{0}_{1}.tif".format(desc1.basename, desc2.basename)
            print nameRaster
            raster_s = (arcpy.Raster(raster_j) - arcpy.Raster(raster_k) )/ (j - k)
            slope1.append(raster_s)
    print "求中值"
    slope_raster = arcpy.sa.CellStatistics(slope1, "MEDIAN", "DATA")
    slope_raster.save(slope_raster_path)

if __name__ == '__main__':
    source_folder = "G:/GPP_Temp/res_mean"
    outpt_folder = "G:/GPP_Temp/MK_res"
    mkSen_cal(source_folder, outpt_folder)

上面的代码中的MK算法与上一篇博客是一样的，总体来看只是加了一个SEN算法。SEN算法中的求中位数我是用了arcpy.sa.CellStatistics(slope1, "MEDIAN", "DATA")来求得，但我在用的时候发现当影像数超过351张后所得到的函数返回影像值均为NODATA。所以又写了下面的代码，下面的代码并不都是程序完成的，需要再借助ARCGIS。

def mkSen_cal(source_folder,outpt_folder,out_slopes):
    arcpy.CheckOutExtension("spatial")
    arcpy.env.workspace = source_folder
    print "源数据的路径为："+str(arcpy.env.workspace)
    arcpy.overwriteOutput = True
    arcpy.CheckOutExtension("spatial")
    muti_path = outpt_folder
    print "成果保存的路径为:"+muti_path
    raster_list = arcpy.ListRasters("*","tif")
    print "tif文件的个数："+str(len(raster_list))
    print raster_list
    sum_array = 0
    n = len(raster_list)

    print "开始求统计量S"
    for j in range(n-1):
        for k in range(j+1,n,1):
        #print str(k)+"--"+str(j)
            raster1 = arcpy.env.workspace + "/"+raster_list[j]
            raster2 = arcpy.env.workspace + "/"+raster_list[k]
            desc1 = arcpy.Describe(raster1)
            desc2 = arcpy.Describe(raster2)
            nameRaster = "正在叠加:"+"{0}_{1}.tif".format(desc1.basename,desc2.basename)
            print nameRaster

            arr = Muti_2page.muti_tif(raster1,raster2)
            sum_array = sum_array + arr
    arcpy.env.workspace = muti_path
    s_Raster_path = "s_Raster.tif"
    print s_Raster_path+"统计量计算成功"
    if arcpy.Exists(s_Raster_path):
        arcpy.Delete_management(s_Raster_path)
        print "正在删除文件......"
    sum_array.save(s_Raster_path)
    print "S统计量栅格计算并保存成功"

    print "开始计算Z统计量"
    s_Raster = arcpy.Raster(muti_path+"/"+"s_Raster.tif")

    var = math.sqrt(n*(n-1)*(2*n+5)/18)
    array1 = arcpy.sa.Con(s_Raster,(s_Raster+1)/var,s_Raster,"VALUE < 0")
    array2 = arcpy.sa.Con(array1,(s_Raster-1)/var,array1,"VALUE > 0")
    z_Raster_path = muti_path+"/"+"z_Raster.tif"
    arcpy.Delete_management(array1)
    if arcpy.Exists(z_Raster_path):
        arcpy.Delete_management(z_Raster_path)
    array2.save(z_Raster_path)
    print "计算Z统计量成功"
    print "开始对Z统计量进行评估"
    Muti_2page.z_belive(z_Raster_path, muti_path)
    print "Z统计量的评估完成"

    arcpy.env.workspace = source_folder
    arcpy.overwriteOutput = True
    arcpy.CheckOutExtension("spatial")
    print "Sen计算倾斜度"
    muti_path = out_slopes
    raster_list = arcpy.ListRasters("*", "tif")
    print "tif文件的个数：" + str(len(raster_list))
    print raster_list
    n = len(raster_list)

    print  "n为：" + str(n)
    for k in range(0, n - 1):
        for j in range(k + 1, n):
            raster_j = raster_list[j]
            raster_k = raster_list[k]
            desc1 = arcpy.Describe(raster_j)
            desc2 = arcpy.Describe(raster_k)
            nameRaster = "正在叠加:" + "{0}_{1}.tif".format(desc1.basename, desc2.basename)
            print nameRaster
            raster_s = (arcpy.Raster(raster_j) - arcpy.Raster(raster_k)) / (j - k)
            path = muti_path + "/" + str(j) + str(k) + ".tif";
            raster_s.save(path)

if __name__ == '__main__':
    source_folder = "G:/GPP_Temp/res_mean"
    outpt_folder = "G:/GPP_Temp/MK_res"
    out_slopes = "G:/GPP_Temp/MK_res/slope"
    mkSen_cal(source_folder, outpt_folder,out_slopes)

这种方法是不用arcpy.sa.CellStatistics（）这个函数，将这一步回到ARCGIS中做，输入数据为G:/GPP_Temp/MK_res/slope"目录下的栅格影像，得成的结果就为SEN结果。

总结：上面的代码得到MK+SEN的两张影像后就可以做接下来的分析了。