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处理气象数据——目录