基于Matlab的栅格水平复相关,偏相关和相关
2020-02-17 本文已影响0人
画长空_yin
在进行影响因素分析时,常使用相关分析,偏相关分析和复相关分析法来量化一种或多种要素对一种因素的影响程度,在前期的博客中分享了关于相关分析和偏相关分析的代码,在本文中我们把三种相关分析集成在一起,同时输出三种结果,该分析是基于长时间的栅格数据进行的,以植被覆盖度,降水和温度为例,具体如下
@ author yinlichang3064@163.com
[a,R]=geotiffread('F:\校级课题项目\data\屏障带\2002降水.tif');%先导入投影信息
info=geotiffinfo('F:\校级课题项目\data\屏障带\2002降水.tif');
[m,n]=size(a);
begin_year=2000; % 开始年份
end_year=2018; % 结束年份
cd=end_year-begin_year+1;
presum=zeros(m*n,cd);
kk=1;
for year=begin_year:end_year
filename=strcat('F:\校级课题项目\data\屏障带\',int2str(year),'降水.tif');
data=importdata(filename);
data=reshape(data,m*n,1);
presum(:,kk)=data;
kk=kk+1;
end
tempsum=zeros(m*n,cd);
kk=1;
for year=begin_year:end_year
filename=strcat('F:\校级课题项目\data\屏障带\',int2str(year),'温度.tif');
data=importdata(filename);
data=reshape(data,size(a,1)*size(a,2),1);
tempsum(:,kk)=data;
kk=kk+1;
end
fcsum=zeros(m*n,cd);
kk=1;
for year=begin_year:end_year
filename=strcat('F:\校级课题项目\data\屏障带\',int2str(year),'植被覆盖度.tif');
data=importdata(filename);
data=reshape(data,m*n,1);
fcsum(:,kk)=data;
kk=kk+1;
end
pc_pre_fvc=zeros(m,n)+nan;
pc_pre_fvc_p=zeros(m,n)+nan;
pc_temp_fvc=zeros(m,n)+nan;
pc_temp_fvc_p=zeros(m,n)+nan;
mpc=zeros(m,n)+nan;
mpc_p=zeros(m,n)+nan;
c_temp_fvc=zeros(m,n)+nan;
c_temp_fvc_p=zeros(m,n)+nan;
c_pre_fvc=zeros(m,n)+nan;
c_pre_fvc_p=zeros(m,n)+nan;
for i=1:m*n
predata=presum(i,:)';
fvcdata=fcsum(i,:)';
tempdata=tempsum(i,:)';
if min(predata)>0 && min(fvcdata)>=0 && min(tempdata)>-100 % 有效值判定0,0,-100
% 相关
[r1,p1]=corrcoef(predata,fvcdata); % 降水和ndvi的相关
c_pre_fvc(i)=r1(2);
c_pre_fvc_p(i)=p1(2);
[r11,p11]=corrcoef(tempdata,fvcdata); % 温度和ndv的相关
c_temp_fvc(i)=r11(2);
c_temp_fvc_p(i)=p11(2);
% 偏相关
[rho,p]=partialcorr(predata,fvcdata,tempdata);%注意,控制的变量放在最后面 控制温度,看降水和植被覆盖度的的偏相关
pc_pre_fvc(i)=rho;
pc_pre_fvc_p(i)=p;
[rho1,p1]=partialcorr(tempdata,fvcdata,predata);%注意,控制的变量放在最后面
pc_temp_fvc(i)=rho1;
pc_temp_fvc_p(i)=p1;
% 复相关 指度量?y \ y?y 与其最优线性预测 y? 的相关系数
X=[ones(cd,1),predata,tempdata];
[b,~,~,~,~] = regress(fvcdata,X);
fvc_yz=b(1)+b(2).*predata+b(3).*tempdata;
[r2,p2]=corrcoef(fvc_yz,fvcdata);
mpc(i)=r2(2);
mpc_p(i)=p2(2);
end
end
result_mulu='D:\';
% 相关系数输出
filename1=[result_mulu,'降水和植被覆盖度相关系数.tif'];
geotiffwrite(filename1,c_pre_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename2=[result_mulu,'降水和植被覆盖度相关系数显著性p值.tif'];
geotiffwrite(filename2,c_pre_fvc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename3=[result_mulu,'通过显著性0.05检验的降水和植被覆盖度相关系数.tif'];
c_pre_fvc(c_pre_fvc_p>0.05)=NaN;
geotiffwrite(filename3,c_pre_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename1=[result_mulu,'温度和植被覆盖度相关系数.tif'];
geotiffwrite(filename1,c_temp_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename2=[result_mulu,'温度和植被覆盖度相关系数显著性p值.tif'];
geotiffwrite(filename2,c_temp_fvc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename3=[result_mulu,'通过显著性0.05检验的温度和植被覆盖度相关系数.tif'];
c_temp_fvc(c_temp_fvc_p>0.05)=NaN;
geotiffwrite(filename3,c_temp_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
% 偏相关系数输出
filename1=[result_mulu,'降水和植被覆盖度的偏相关系数_控制温度.tif'];
geotiffwrite(filename1,pc_pre_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename2=[result_mulu,'降水和植被覆盖度的偏相关系数显著性p值_控制温度.tif'];
geotiffwrite(filename2,pc_pre_fvc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename3=[result_mulu,'通过显著性0.05检验的降水和植被覆盖度的偏相关系数_控制温度.tif'];
pc_pre_fvc(pc_pre_fvc_p>0.05)=NaN;
geotiffwrite(filename3,pc_pre_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename4=[result_mulu,'温度和植被覆盖度的偏相关系数_控制降水.tif'];
geotiffwrite(filename4,pc_temp_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename5=[result_mulu,'温度和植被覆盖度的偏相关系数显著性p值_控制降水.tif'];
geotiffwrite(filename5,pc_temp_fvc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename6=[result_mulu,'通过显著性0.05检验的温度和植被覆盖度的偏相关系数_控制降水.tif'];
pc_temp_fvc(pc_temp_fvc_p>0.05)=NaN;
geotiffwrite(filename6,pc_temp_fvc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
% 复相关系数输出
filename7=[result_mulu,'温度降水与植被覆盖度的复相关系数.tif'];
geotiffwrite(filename7,mpc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename8=[result_mulu,'温度降水与植被覆盖度的复相关系数显著性p值.tif'];
geotiffwrite(filename8,mpc_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
filename9=[result_mulu,'通过显著性0.05检验的温度降水与植被覆盖度的复相关系数.tif'];
mpc(mpc_p>0.05)=NaN;
geotiffwrite(filename9,mpc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
上述代码会输出相关系数,相关系数的显著性p值(F检验),以及通过显著性检验的相关系数,如有更多需求可邮件联系