MK检验_MK趋势检验_sen+mk_

clear [a,R]=geotiffread('D:\GIS\vegetation\output\yearmax\1982.tif'); %先导入投影信息 info=geotiffinfo('D:\GIS\vegetation\output\yearmax\1982.tif');%先导入投影信息 [m,n]=size(a); cd=34; %34年，时间跨度 datasum=zeros(m*n,cd)+NaN; p=1; for year=1982:2015 %起始年份 filename=['D:\qixiang\年全国8kmPET\china',int2str(year),'pet.tif']; data=importdata(filename); data=reshape(data,m*n,1); datasum(:,p)=data; % p=p+1; end invalid_value=datasum(1); sresult=zeros(m,n)+NaN; for i=1:size(datasum,1) % data=datasum(i,:); if max(data)>invalid_value % sgnsum=[]; for k=2:cd for j=1:(k-1) sgn=data(k)-data(j); if sgn>0 sgn=1; else if sgn<0 sgn=-1; else sgn=0; end end sgnsum=[sgnsum;sgn]; end end add=sum(sgnsum); sresult(i)=add; end end vars=cd*(cd-1)*(2*cd+5)/18; zc=zeros(m,n)+NaN; sy=find(sresult==0); zc(sy)=0; sy=find(sresult>0); zc(sy)=(sresult(sy)-1)./sqrt(vars); sy=find(sresult<0); zc(sy)=(sresult(sy)+1)./sqrt(vars); geotiffwrite('C:\MATLAB\MK检验结果.tif',zc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag); %注意修改路径 通过上述代码的运行可以得到MK检验的结果。上述代码运行时只需要修改起始年份和年份长度以及文件的名称，注意文件名称 按照规律来进行分布，本文中的名称是china1982pet.tif，china1983pet.tif...china2015pet.tif,保证能够按照规律读取。 假设读者已经运行完了sen代码和本文中的代码，则可以得到两个tif文件，分别是MK检验结果和sen的结果，进而通过以下代码 来进行最终的判断 [a,R]=geotiffread('D:\GIS\vegetation\output\yearmax\1982.tif'); %先导入投影信息 info=geotiffinfo('D:\GIS\vegetation\output\yearmax\1982.tif');%先导入投影信息 data=importdata('C:\MATLAB\MK检验结果.tif'); sen_value=importdata('D:\zhang\基于sen的pet变化趋势.tif'); sen_value(abs(data)<1.94)=NaN; %MK结果值高于1.94则认为通过了显著性95% geotiffwrite('C:\MATLAB\通过显著性95%的MK+sen趋势分析结果.tif',sen_value,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);