Sen+MK长时间序列趋势性分析----基于python的代码实现

最新推荐文章于 2024-08-21 16:26:07 发布

P0级推荐算法工程师

最新推荐文章于 2024-08-21 16:26:07 发布

阅读量1.4w

点赞数 75

分类专栏： python空间数据处理文章标签： python numpy

本文链接：https://blog.csdn.net/sunstronginchina/article/details/115624002

版权

python空间数据处理专栏收录该内容

8 篇文章 16 订阅

订阅专栏

sen+mk python实现代码免费共享-----赶紧收藏吧

python开源社区公布了进行sen+mk趋势性检验的官方包，有关该官方包的主要内容详见：https://github.com/Coder2cdb/pyMannKendall，开源不易，有些人还收费查看，我觉得没必要，本次代码免费分享给大家，有需要的收藏点赞！

本次sen+mk python代码输入参数：多年影像的存储路径
path1=r"E:\2021-03-03—生态修复\吴起本地\NDVI_wuqi"
本次sen+mk python代码输出参数：结果影像存储路径
result_path=r"E:\生态修复\1990-2000_NDVI\result"
调用方法：
sen_mk_test(path1, result_path)
MannKendall输出参数说明

''       
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
'''

1.导入python包

#coding:utf-8
import numpy as np
import pymannkendall as mk
import os 
import rasterio as ras

2.获取所有的影像路径，按照从小到大年份排序

path1=r"E:\2021-03-03—生态修复\吴起本地\NDVI_wuqi"
filepaths=[]
for file in os.listdir(path1):
    filepath1=os.path.join(path1,file)
    filepaths.append(filepath1)

3.读取所有的影像数据并拼接为一个numpy矩阵

#获取影像数量
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

4.定义输出矩阵

#输出矩阵，无值区用-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检验，如果所有影像同一像元都为空，则不进行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])))

5.sen+mk检验

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)
    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

6.定义写影像的函数

#写影像，包括宽度，高度，投影，波段名，矩阵
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

7.输出结果

#输出矩阵
all_array=[slope_array,Trend_array,p_array,s_array,Tau_array,z_array]   
#输出路径
result_path=r"E:\2021-03-03—生态修复\吴起本地\result"
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)

sen+mk代码汇总

'''
Created on 2021年4月11日

@author: SunStrong
'''
#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:结果输出路径
    
    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:\2021-03-03—生态修复\吴起本地\NDVI_wuqi"
result_path=r"E:\2021-03-03—生态修复\吴起本地\result"
sen_mk_test(path1, result_path)