MATLAB--按显著性分级的相关系数空间分布图
作者：第八星系-工程师
邮箱：2921725518@qq.com
数据处理
%以青藏高原81-97年地表径流与2m气温为例
ID1 = ['F:\毕设\era5\81-19\data.nc']; %数据路径
sro1 = ncread(ID1,'sro'); %读取变量地表径流
t2m1 = ncread(ID1,'t2m'); %读取变量2m气温
lat = ncread(ID1,'latitude'); %读取纬度信息
lon = ncread(ID1,'longitude'); %读取经度信息
%将月数据以年为单位排放
k = 1;
for i = 1:39
msro(:,:,:,i) = sro1(:,:,k:k+11);
mt2m(:,:,:,i) = t2m1(:,:,k:k+11);
k = k+12;
end
%求取81-97年年平均
ysro = squeeze(mean(msro(:,:,:,1:17),3));
yt2m = squeeze(mean(mt2m(:,:,:,1:17),3));计算相关系数
for i = 1:301
for j = 1:151
t2m=yt2m(i,j,:);
t2m = squeeze(t2m); %将数组压缩为17*1
sro=ysro(i,j,:);
sro = squeeze(sro);
r2 = corr(t2m,sro); %用corr函数求两个变量的相关系数
r(i,j) = r2;
end
end固定目标区域
%读取高原边界经纬度数据
Boundary = load('F:\剑宗密史\冻融指数\指数值340x180x49\TP_boundary.txt');
V1 = Boundary(:,1);
V2 = Boundary(:,2);
%屏蔽高原以外的数据
for i = 1:301
for j = 1:151
in = inpolygon(lon(i),lat(j),V1,V2);
if in == 0
r(i,j) = NaN;
end
end
end
r1 = r;显著性等级分级
%相关系数显著性检验表
%https://www.docin.com/p-272983879.html
%此处1为显著性达95%的负相关；2为显著性达90%的负相关；3为显著性90%以下负相关；4为显著性90%以下正相关；5为显著性达90%的正相关；6为显著性达95%的正相关
%17 0.389 0.456
for i = 1:301
for j = 1:151
if r(i,j)>0.456 %令相关系数大于0.456的位置都等于6（达95%显著性）
r1(i,j) = 6;
elseif r(i,j)>0.389&&r(i,j)<0.456 %令相关系数位于0.389和0.456之间位置都等于5（达90%显著性）
r1(i,j) = 5;
elseif r(i,j)<-0.456
r1(i,j) = 1;
elseif r(i,j)>-0.456&&r(i,j)<-0.389
r1(i,j) = 2;
elseif 0<r(i,j)&&r(i,j)<0.389
r1(i,j) = 4;
elseif -0.389<r(i,j)&&r(i,j)<0
r1(i,j) = 3;
end
end
end
%% 统计各显著性的频数和频率
P1 = 0;
P2 = 0;
P3 = 0;
N1 = 0;
N2 = 0;
N3 = 0;
%令各计数变量满足某显著性的位置等于1
for i = 1:301
for j = 1:151
if r1(i,j) == 4 %如果(i,j)位置上r1等于4（显著性90%以下的正相关）P1计数+1，其他显著等级类似
P1 = P1 + 1;
elseif r1(i,j) == 5
P2 = P2 + 1;
elseif r1(i,j) == 6
P3 = P3 + 1;
elseif r1(i,j) == 1
N3 = N3 + 1;
elseif r1(i,j) == 2
N2 = N2 + 1;
elseif r1(i,j) == 3
N1 = N1 + 1;
end
end
end
%所得P1、P2、P3、N1、N2、N3即为各显著等级频数
ngrad = P1+P2+P3+N1+N2+N3; %总有效格点数（没有r=0的情况下）
%各显著性的频率
f(1,1) = P1/ngrad;
f(2,1) = P2/ngrad;
f(3,1) = P3/ngrad;
f(4,1) = N3/ngrad;
f(5,1) = N2/ngrad;
f(6,1) = N1/ngrad;作图
figure(1)
contourf(lon,lat,r1','LineStyle','none'); %用contourf函数绘制等值线图，此处设置了关闭等值线
hold on
axis() %确定x轴与y轴框图大小
set(gca,'YTick',) %y轴刻度
set(gca,'XTick',) %x轴刻度
set(gca,'yticklabel',{'25°N','30°N','35°N','40°N'})%y轴刻度标签
set(gca,'xticklabel',{'75°E','80°E','85°E','90°E','90°E','95°E','100°E','105°E'})%x轴刻度标签
set(gca,'CLim',) %等值线数值范围
hold on
map = ;
colormap(map) %map是自定义rgb色卡，用于等值线图填色
title('1981-1997年青藏高原地表径流与2m气温相关系数') %标题
hold on
plot(V1,V2,'k','LineWidth',1) %高原边界
colorbar('Position',,'Ticks',[],'TickLabels',{''});%设置色标位置，关闭刻度和刻度标签
%在色标旁边标注对应显著程度，实际位置可在图窗调整
annotation(figure(1),'textbox',,'String',{'P**'},'FitBoxToText','off','EdgeColor','none');
annotation(figure(1),'textbox',,'String',{'P*'},'FitBoxToText','off','EdgeColor','none');
annotation(figure(1),'textbox',,'String',{'P'},'FitBoxToText','off','EdgeColor','none');
annotation(figure(1),'textbox',,'String',{'N'},'FitBoxToText','off','EdgeColor','none');
annotation(figure(1),'textbox',,'String',{'N*'},'FitBoxToText','off','EdgeColor','none');
annotation(figure(1),'textbox',,'String',{'N**'},'FitBoxToText','off','EdgeColor','none');
%显著性频率直方图
hold on
x = 1:6;
axes1 = axes('Position',);%设置直方图位置和大小
hold(axes1,'on');
%每种显著性用和等值线图对应的颜色
b1 = bar(1,f(1,1),'FaceColor',map(4,:),'EdgeColor','none','BarWidth',0.4);
b2 = bar(2,f(2,1),'FaceColor',map(5,:),'EdgeColor','none','BarWidth',0.4);
b3 = bar(3,f(3,1),'FaceColor',map(6,:),'EdgeColor','none','BarWidth',0.4);
b4 = bar(4,f(4,1),'FaceColor',map(1,:),'EdgeColor','none','BarWidth',0.4);
b5 = bar(5,f(5,1),'FaceColor',map(2,:),'EdgeColor','none','BarWidth',0.4);
b6 = bar(6,f(6,1),'FaceColor',map(3,:),'EdgeColor','none','BarWidth',0.4);
box(axes1,'off');
set(axes1,'XTick',,'XTickLabel',{'','','','','',''});%直方图刻度及其标签
微信搜索“第八星系人造大气理论爱好者”公众号，关注获取文章数据
页:
[1]