基于matlab的栅格数据相关分析及显著性检验

相关分析在科研过程中常常会碰到，表征了一种两个因素间的相关程度，其值的大小只能说明相关性的强弱。
比如当我们有2000-2015年的NPP数据和2000-2015年的降水数据时，我们想查看两者在空间上随时间变化的
相关性。本文以三者间的相关性为例，分别是产水，NPP和土壤保持
具体代码如下：
%authour@yinlichang3064@163.com

%将三者多年的数据放在三个不同的矩阵中
nppsum=zeros(3587*4642,16);
for year=2000:2015
    filename=strcat('F:\课题项目\data\',int2str(year),'npp.tif');
    data=importdata(filename);
    data=reshape(data,3587*4642,1);
    nppsum(:,year-1999)=data;
end

scsum=zeros(3587*4642,16);
for year=2000:2015
    filename=strcat('F:\课题项目\data\',int2str(year),'sc.tif');
    data=importdata(filename);
    data=reshape(data,3587*4642,1);
    scsum(:,year-1999)=data;
end

wcsum=zeros(3587*4642,16);
for year=2000:2015
    filename=strcat('F:\课题项目\data\',int2str(year),'water_yield.tif');
    data=importdata(filename);
    data=reshape(data,3587*4642,1);
    wcsum(:,year-1999)=data;
end
%相关性和显著性
sc_npp_xgx=zeros(3587,4642);
sc_npp_p=zeros(3587,4642);
sc_wc_xgx=zeros(3587,4642);
sc_wc_p=zeros(3587,4642);
npp_wc_xgx=zeros(3587,4642);
npp_wc_p=zeros(3587,4642);
for i=1:length(scsum)
    sc=scsum(i,:);
    if min(sc)>=0 %注意这里的土壤保持的有效范围是大于0，如果自己的数据有效范围有小于0的话，则可以不用加这个
        %npp与sc
        npp=nppsum(i,:);
        [r,p]=corrcoef(sc,npp);
         sc_npp_xgx(i)=r(2);
         sc_npp_p(i)=p(2);
         wc=wcsum(i,:);
         [r1,p1]=corrcoef(sc,wc);
         sc_wc_xgx(i)=r1(2);
         sc_wc_p(i)=p1(2);
         [r2,p2]=corrcoef(npp,wc);
         npp_wc_xgx(i)=r2(2);
         npp_wc_p(i)=p2(2);
    end
end
filename1='F:\result\sc_npp相关性.tif';
filename2='F:\result\sc_npp显著性.tif';
filename3='F:\sc_wc相关性.tif';
filename4='F:\result\sc_wc显著性.tif';
filename5='F:\result\npp_wc相关性.tif';
filename6='F:\npp_wc显著性.tif';

[a,R]=geotiffread('F:\项目\data\\2002water_yield.tif');%先导入投影信息
info=geotiffinfo('F:\项目\data\2002water_yield.tif');
%%输出图像
geotiffwrite(filename1,sc_npp_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename2,sc_npp_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename3,sc_wc_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename4,sc_wc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename5,npp_wc_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename6,npp_wc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

如果自己的数据中存在负值，不想加限制条件的话，可以用下面的代码，但运行效率会降低，代码如下

for i=1:length(scsum)
    sc=scsum(i,:);
        %npp与sc
        npp=nppsum(i,:);
        [r,p]=corrcoef(sc,npp);
         sc_npp_xgx(i)=r(2);
         sc_npp_p(i)=p(2);
         wc=wcsum(i,:);
         [r1,p1]=corrcoef(sc,wc);
         sc_wc_xgx(i)=r1(2);
         sc_wc_p(i)=p1(2);
         [r2,p2]=corrcoef(npp,wc);
         npp_wc_xgx(i)=r2(2);
         npp_wc_p(i)=p2(2);
end
filename1='F:\result\sc_npp相关性.tif';
filename2='F:\result\sc_npp显著性.tif';
filename3='F:\sc_wc相关性.tif';
filename4='F:\result\sc_wc显著性.tif';
filename5='F:\result\npp_wc相关性.tif';
filename6='F:\npp_wc显著性.tif';

[a,R]=geotiffread('F:\项目\data\\2002water_yield.tif');%先导入投影信息
info=geotiffinfo('F:\项目\data\2002water_yield.tif');
%%输出图像
geotiffwrite(filename1,sc_npp_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename2,sc_npp_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename3,sc_wc_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename4,sc_wc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename5,npp_wc_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename6,npp_wc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);