sen+mk python实现代码免费共享-----赶紧收藏吧
python开源社区公布了进行sen+mk趋势性检验的官方包，有关该官方包的主要内容详见：/Coder2cdb/pyMannKendall，开源不易，有些人还收费查看，我觉得没必要，本次代码免费分享给大家，有需要的收藏点赞！
本次sen+mk python代码输入参数：多年影像的存储路径
path1=r"E:\-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 testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: Theil-Sen estimator/slopeintercept: intercept of Kendall-Theil Robust Line'''
1.导入python包
#coding:utf-8import numpy as npimport pymannkendall as mkimport os import rasterio as ras
2.获取所有的影像路径，按照从小到大年份排序
path1=r"E:\-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.transformheight1=img1.heightwidth1=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=-1elif trend=="increasing":trend_value=1else:trend_value=0slope_array[x,y]=slope#senslopes_array[x,y]=sz_array[x,y]=zTrend_array[x,y]=trend_valuep_array[x,y]=pTau_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:\-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 4月11日@author: SunStrong'''#coding:utf-8import numpy as npimport pymannkendall as mkimport os import rasterio as rasdef 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.transformheight1=img1.heightwidth1=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 testfor 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 testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: Theil-Sen estimator/slopeintercept: intercept of Kendall-Theil Robust Line'''if trend=="decreasing":trend_value=-1elif trend=="increasing":trend_value=1else:trend_value=0slope_array[x,y]=slope#senslopes_array[x,y]=sz_array[x,y]=zTrend_array[x,y]=trend_valuep_array[x,y]=pTau_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:\-03-03—生态修复\吴起本地\NDVI_wuqi"result_path=r"E:\-03-03—生态修复\吴起本地\result"sen_mk_test(path1, result_path)
