MATLAB两类栅格数据之间的相关性计算

两类栅格数据之间的相关性计算（输出为tif影像）

栅格数据做相关性分析前的预处理（批量定义投影、栅格投影、重采样）

栅格影像行列号需要一致，行列号不一致可以在ArcGIS中批量处理：
1.重采样
2.裁剪 右键空白 环境设置


处理完毕后，进行相关分析。

Matlab代码-年

[a,R]=geotiffread('F:\SMC_V3.0\std_anlisyt\annual\x\FP_FVC_2013.8_class.tif_clp.tif');%先导入投影信息，某个影像的路径就行（最好是你分析的数据中的一个）
info=geotiffinfo('F:\SMC_V3.0\std_anlisyt\annual\x\FP_FVC_2013.8_class.tif_clp.tif');%同上
[m,n]=size(a);

nppsum=zeros(m*n,8);%此处要修改，共几年就填写多少，我这里是8年的
for year=2013:2020
    filename=strcat('F:\yz\FVC-time\2013-2020-FVC-time\Re-FVC-time\FVC_farm_perfect-clip\','FP_FVC_',int2str(year),'.8_class.tif_clp.tif');%此处要修改，我这里是八年的每年年均植被覆盖度的数据，注意你的文件名字。
    data=importdata(filename);
    data=reshape(data,m*n,1);
    nppsum(:,year-2012)=data;%此处需要修改，我的数据是从2013开始，此处就为2012.
end

wcsum=zeros(m*n,8);
for year=2013:2020
    filename=strcat('F:\yz\YZ-Date\temperature_date\MOD11B3\autumn(9,10,11)\sample_clip\',int2str(year),'244.sample.tif_clp.tif');%此处要修改，我这里是八年的每年年均地表温度的数据，注意你的文件名字。
    data=importdata(filename);
    data=reshape(data,m*n,1);
    wcsum(:,year-2012)=data;
end

%相关性和显著性
npp_wc_xgx=zeros(m,n);
npp_wc_p=zeros(m,n);
for i=1:length(nppsum)
    npp=nppsum(i,:);
    if min(npp)>0 %注意这里的NPP的有效范围是大于0，如果自己的数据有效范围有小于0的话，则可以不用加这个
        wc=wcsum(i,:);
         [r2,p2]=corrcoef(npp,wc);
         npp_wc_xgx(i)=r2(2);
         npp_wc_p(i)=p2(2);
    end
end

filename5='F:\yz\YZ-Date\temperature_date\MOD11B3\autumn(9,10,11)\std\相关性_LST秋季节尺度.tif';%此处要修改，输出的路径及名字
filename6='F:\yz\YZ-Date\temperature_date\MOD11B3\autumn(9,10,11)\std\显著性_LST秋季节尺度.tif';%同上


%%输出图像
geotiffwrite(filename5,npp_wc_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename6,npp_wc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

参考程序来源于：参考链接

Matlab代码-月

%先导入投影信息，某个影像的路径就行（最好是你分析的数据中的一个）
[a,R]=readgeoraster('G:\SIF\Global-AI_monthly_v3\bi\199001.tif');
info=geotiffinfo('G:\SIF\Global-AI_monthly_v3\bi\199001.tif');
[m,n]=size(a);
i=1;
 
gw=zeros(m*n,24); 
%此处要修改，共多少个月就写多少个月,我这里一共51个月，gw是我自己的变量名，可以根据需求修改
 
    for year=1990:1991
        for month=1:12
        filename=strcat('G:\SIF\Global-AI_monthly_v3\bi\',int2str(year),num2str(month,'%02d'),'.tif');
        %上述是文件名字的写法，我这里的例子是YYYYMM，根据自己文件命名修改
        data=importdata(filename);
        data=reshape(data,m*n,1);
        gw(:,i)=data; 
        i=i+1;
        end
    end
   
    % 可以看出每一年只挑选3-5月份的数据，用的双层循环，如果是每年的十二个月，就month=1:12
 
i=1;  %前面一个参数，下面导入另一个参数rssm，这里需要重新命i=1!
rssm=zeros(m*n,24);
   for year=1990:1991
        for month=1:12
        filename=strcat('G:\GPP_he\npp\test\',int2str(year),num2str(month,'%02d'),'.tif');
        data=importdata(filename);
        data=reshape(data,m*n,1);
        rssm(:,i)=data; 
        i=i+1;
        end
   end
%求相关性和显著性
 
grw_sm_xgx=zeros(m,n);
grw_sm_p=zeros(m,n);
for tl=1:length(gw)     %tl=timeline 时间长度，其实就是前面的51.
    grw=gw(tl,:);
    ssm=rssm(tl,:);
    [r2,p2]=corr(grw',ssm','type',"Pearson"); 
    %这里可以选自己想要的相关性，有三种。常用Pearson
    grw_sm_xgx(tl)=r2;
    grw_sm_p(tl)=p2;
end
 
filename1='G:\SIF\Global-AI_monthly_v3\result\相关性.tif';   %相关性
filename2='G:\SIF\Global-AI_monthly_v3\result\显著性.tif';   %显著性
 
%输出图像
geotiffwrite(filename1,grw_sm_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename2,grw_sm_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

参考