Matlab处理气象数据（十五）分省份的趋势分析

用省份边界数据，制作掩膜（mask），从而提取出各省的DTR数据。并且求P和R2。

%求P值
t=1:35;
t=t';
t=[ones(35,1),t];
[b,bint,r,rint,stats]= regress(NCEP,t,0.05);

%结果：
alpha为显著性水平(缺省时设定为0.05)

stats是用于检验回归模型的统计量有四个数值
R2，其中R是相关系数
F统计量值
概率P，与统计量F对应
an estimate of the error variance（一个错误的方差估计）

%回归系数估计值
b =

-0.5499
0.0306

%b的置信区间
bint =

-0.8822 -0.2177
0.0145 0.0466

%残差
r =

0.7600
-0.3664
-0.2539
0.0911
-0.0848
-0.6574
-0.3917
-0.3051
-0.0591
-0.1611
-0.2956
0.2641
-0.4424
-0.0424
-0.2965
0.8728
-0.0038
-0.1784
0.4066
0.6409
0.1986
0.3052
0.4790
0.8306
-0.1097
0.5178
-0.0808
0.5138
0.7762
-0.5602
-0.0756
-0.4685
-0.8054
-0.7835
-0.2346

%r的置信区间
rint =

-0.1201 1.6401
-1.2834 0.5505
-1.1802 0.6724
-0.8433 1.0256
-1.0232 0.8535
-1.5696 0.2549
-1.3268 0.5433
-1.2473 0.6372
-1.0101 0.8919
-1.1130 0.7908
-1.2455 0.6542
-0.6887 1.2170
-1.3882 0.5034
-1.0026 0.9179
-1.2519 0.6589
-0.0366 1.7821
-0.9662 0.9587
-1.1388 0.7821
-0.5446 1.3579
-0.2931 1.5749
-0.7601 1.1573
-0.6488 1.2593
-0.4645 1.4225
-0.0792 1.7403
-1.0647 0.8452
-0.4175 1.4531
-1.0317 0.8700
-0.4166 1.4443
-0.1273 1.6797
-1.4808 0.3604
-1.0140 0.8629
-1.3882 0.4512
-1.6900 0.0792
-1.6660 0.0990
-1.1523 0.6831

%stats是用于检验回归模型的统计量有四个数值
第一个是R2，其中R是相关系数，第二个是F统计量值，第三个是与统计量F对应的概率P，第四个是 an estimate of the error variance（一个错误的方差估计）。
stats =

0.3113 14.9132 0.0005 0.2234

说明：
R2表示方差解释率，R2越接近1说明数据拟合程度越好。
F统计量用于检验模型是否通过检验。通过查F分布表，如果F>F分布表中对应的值，则通过检验。

P为F 统计量对应的概率，越接近0越好，当P<α时拒绝H0，回归模型成立。

以下程序以安徽为例，其他省份只需将“安徽”替换为其他省份的名称即可。

%% 
%抠NCEP数据
clc;clear
lon=72:0.5:135.5;
lat=18:0.5:53.5;
[x y]=meshgrid(lon,lat); 
 
load('E:\数据\201501\T1')
T10=reshape(T1,[72 128 12 35]);
T11=reshape(T1,[72 128 12*35]);
T10=squeeze(nanmean(T10,3));
T_mean=squeeze(nanmean(T11,3));
clear T1
for i=1:35
T1(:,:,i)=squeeze(T10(:,:,i))-T_mean;
end
map=shaperead('E:\数据\边界\各省区划边界\安徽\安徽.shp');
bou1_4lx=[map(:).X];%提取经度信息 
bou1_4ly= [map(:).Y];%提取纬度信息 
R=makerefmat('RasterSize',size(T1),'Lonlim',[72 135.5],'Latlim',[18 53.5]);
MASK=vec2mtx(bou1_4ly,bou1_4lx,squeeze(T1(:,:,1)),R,'filled');
MASK(find(MASK>1))=nan;
MASK(find(MASK==1))=0;
T12=T1;%删除边界区域外的值
for i=1:72
    for j=1:128
if isnan(MASK(i,j))==1
T12(i,j,:)=nan;
end
    end
end
 
% clf
% contourf(x,y,squeeze(T2(:,:,1)),30);
% shading flat
% colorbar
% hold on
% plot(bou1_4lx,bou1_4ly,'-k','linewidth',3)
% hold off
 
T12=reshape(T12,[72*128 35]);
T12=nanmean(T12,1);
NCEP=T12;
p1=polyfit(1:35,T12,1);  % p输出两个值第一个是a  第二个是b   y=ax+b
target1=polyval(squeeze(p1(:)),1:35);
clear y x lon lat map j i bou1_4ly bou1_4lx  T1 MASK R T12 T10 T11 T_mean
 
%%
%抠观测数据
lon=72:0.5:135.5;
lat=18:0.5:53.5;
[x y]=meshgrid(lon,lat); 
% [lon lat]=meshgrid([97:0.1:107],[21:0.1:30]); 
load('E:\数据\201501\T2')
T10=reshape(T2,[72 128 12 35]);
T11=reshape(T2,[72 128 12*35]);
T10=squeeze(nanmean(T10,3));
T_mean=squeeze(nanmean(T11,3));
clear T2
for i=1:35
T2(:,:,i)=squeeze(T10(:,:,i))-T_mean;
end
map=shaperead('E:\数据\边界\各省区划边界\安徽\安徽.shp');
bou1_4lx=[map(:).X];%提取经度信息 
bou1_4ly= [map(:).Y];%提取纬度信息 
R=makerefmat('RasterSize',size(T2),'Lonlim',[72 135.5],'Latlim',[18 53.5]);
MASK=vec2mtx(bou1_4ly,bou1_4lx,squeeze(T2(:,:,1)),R,'filled');
MASK(find(MASK>1))=nan;
MASK(find(MASK==1))=0;
T22=T2;
for i=1:72
    for j=1:128
if isnan(MASK(i,j))==1
T22(i,j,:)=nan;
end
    end
end
 
% clf
% contourf(x,y,squeeze(T3(:,:,1)),30);
% shading flat
% colorbar
% hold on
% plot(bou1_4lx,bou1_4ly,'-k','linewidth',3)
% hold off
 
T22=reshape(T22,[72*128 35]);
T22=nanmean(T22,1);
Obs=T22;
p2=polyfit(1:35,T22,1);  % p输出两个值第一个是a  第二个是b   y=ax+b
target2=polyval(squeeze(p2(:)),1:35);
clear y x lon lat map j i bou1_4ly bou1_4lx T2 T1 MASK R  T22 T10 T11 T_mean
%%
%M-k test
r=zeros(1,35);
for k=1:35
      for l=1:k
        if NCEP(k)>=NCEP(l)
     r(k)= r(k) +1;
        end
     end
end
sk=zeros(1,35);
Ek=zeros(1,35);
var_sk=zeros(1,35);
for k=2:35
sk(k)=sum(r(1:k));
end
for k=1:35;
Ek(k)=k*(k-1)/4;
var_sk(k)=k*(k-1)*(2*k+5)/72;
end
for k=1:35;
UFk1(k)=(sk(k)-Ek(k))/sqrt(var_sk(k));
end
clear Ek UFk k l r sk var_sk
%%%%%%%%%%%%%%%%%%%%%%%%
r=zeros(1,35);
for k=1:35
      for l=1:k
        if Obs(k)>=Obs(l)
     r(k)= r(k) +1;
        end
     end
end
sk=zeros(1,35);
Ek=zeros(1,35);
var_sk=zeros(1,35);
for k=2:35
sk(k)=sum(r(1:k));
end
for k=1:35;
Ek(k)=k*(k-1)/4;
var_sk(k)=k*(k-1)*(2*k+5)/72;
end
for k=1:35;
UFk2(k)=(sk(k)-Ek(k))/sqrt(var_sk(k));
end
clear Ek UFk k l r sk var_sk
 
r0=corrcoef(squeeze(Obs(:)),squeeze(NCEP(:)));%求Obs和NCEP的相关性
R=r0(1,2);
%%
clf %clear figure清理图像
subplot(2,1,1) %图两行一列的第一张图
plot(Obs,'b','linewidth',2)
hold on
plot(NCEP,'r','linewidth',2)
plot(target1,'r','linewidth',2)
plot(target2,'b','linewidth',2)
 
legend({'Observed','NCEP'},'Location','Northwest')
set(gca,'xtick',[2 7 12 17 22 27 32]); 
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'}); 
axis([ -inf inf -1.5 1.5]);
h=xlabel('Year');
set(h, 'FontSize',8,'FontWeight','Bold')
h=ylabel('Temp. Anom. (\circC)');
set(h, 'FontSize',8,'FontWeight','Bold')
set(gca, 'FontSize',8,'FontWeight','Bold','tickdir','out')
xlim([1 35])%x轴范围锁定为1~35
box off %去掉外框
h=text(31,-0.7,['R='  num2str(R) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
aa=zeros(1,35);
bb=aa;
aa(:)=1.64485;
bb(:)=-1.64485;
subplot(2,1,2)%第二张图
plot(UFk2,'b','linewidth',2)
hold on
plot(UFk1,'r','linewidth',2)
plot(aa,'k--','linewidth',2)
plot(bb,'k--','linewidth',2)
legend({'Observed','NCEP'},'Location','Northwest')
set(gca,'xtick',[2 7 12 17 22 27 32]); 
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'}); 
h=xlabel('Year')
set(h, 'FontSize',8,'FontWeight','Bold');
xlim([1 35])
set(gca, 'FontSize',8,'FontWeight','Bold','tickdir','out')
box off

得到的各省的趋势分析如图：

我国31个省、市、自治区的年均DTR异常

---

相关链接：
Matlab处理气象数据——目录