1200字范文,内容丰富有趣,写作的好帮手!
1200字范文 > mk趋势分析matlab 基于matlab 的长时间栅格数据的Sen+MK显著性检验趋势分析

mk趋势分析matlab 基于matlab 的长时间栅格数据的Sen+MK显著性检验趋势分析

时间:2022-12-21 11:32:08

相关推荐

mk趋势分析matlab 基于matlab 的长时间栅格数据的Sen+MK显著性检验趋势分析

在前一篇文章中讲述了用sen法进行长时间的趋势分析,但并未对结果进行显著性检验,通常Sen与MK检验是结合在一起的,

因此本文主要讲述如何进行MK检验。具体代码如下

% @author yinlichang3064@

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: %起始年份

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

sresult=zeros(m,n)+NaN;

for i=1:size(datasum,1) %

data=datasum(i,:);

if min(data)>0 % 有效格点判定,我这里有效值在0以上

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...chinapet.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.96)=NaN; %MK结果值高于1.96则认为通过了显著性95%

geotiffwrite('C:\MATLAB\通过显著性95%的MK+sen趋势分析结果.tif',sen_value,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);%注意修改路径

更多需求,请查看个人介绍

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。