技術標籤：北京二年matlab筆記

在前一篇文章中講述了用sen法進行長時間的趨勢分析，但並未對結果進行顯著性檢驗，通常Sen與MK檢驗是結合在一起的，
因此本文主要講述如何進行MK檢驗。具體程式碼如下

% @author [email protected]
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
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...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.96)=NaN; %MK結果值高於1.96則認為通過了顯著性95%
geotiffwrite('C:\MATLAB\通過顯著性95%的MK+sen趨勢分析結果.tif',sen_value,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);%注意修改路徑