Matlab处理⽓象数据（⼗七）季节变化17.1 各省冬季和夏季的DTR趋势
以江西为例，其他省份只需将下⾯代码中的“江西”替换为其他省份的名称即可。
%%
clc,clear
load('E:\数据\201507\各省MASK\江西');%MASK是各省掩膜
load('E:\数据\201501\Tmax1');%导⼊NCEP数据最⾼温度72*128*396
load('E:\数据\201501\Tmax2');%导⼊观测数据最⾼温度72*128*396
load('E:\数据\201501\Tmin1');%导⼊NCEP数据最低温度72*128*396
load('E:\数据\201501\Tmin2');%导⼊观测数据最低温度72*128*396
DTR1=Tmax1-Tmin1;
DTR2=Tmax2-Tmin2;
%冬季DTR温度72*128*32
for i=1:32
wDTR1(:,:,i)=squeeze(nanmean(DTR1(:,:,(12+12*(i-1)):(14+12*(i-1))),3));
wDTR2(:,:,i)=squeeze(nanmean(DTR2(:,:,(12+12*(i-1)):(14+12*(i-1))),3));
end
%全国夏季DTR温度72*128*33
for i=1:33
smDTR1(:,:,i)=squeeze(nanmean(DTR1(:,:,(6+12*(i-1)):(7+12*(i-1))),3));
smDTR2(:,:,i)=squeeze(nanmean(DTR2(:,:,(6+12*(i-1)):(7+12*(i-1))),3));
end
clear Tmax1 Tmax2 Tmin1 Tmin2 DTR1 DTR2
%%
%冬季
%对NCEP数据
T_mean=squeeze(nanmean(wDTR1,3));
for i=1:32
wDTR1(:,:,i)=squeeze(wDTR1(:,:,i))-T_mean;
end
T12=wDTR1;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T12(i,j,:)=nan;
end
end
end
T12=reshape(T12,[72*12832]);
T12=nanmean(T12,1);
NCEP=T12;
p1=polyfit(1:32,T12,1);% p输出两个值第⼀个是a  第⼆个是b  y=ax+b
target1=polyval(squeeze(p1(:)),1:32);
clear wDTR1 T12 T_mean
%对观测数据
T_mean=squeeze(nanmean(wDTR2,3));
for i=1:32
wDTR2(:,:,i)=squeeze(wDTR2(:,:,i))-T_mean;
end
clear i j T_mean
T22=wDTR2;%删除边界区域外的值
for i=1:72
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T22(i,j,:)=nan;
end
end
end
T22=reshape(T22,[72*12832]);
T22=nanmean(T22,1);
Obs=T22;
p2=polyfit(1:32,T22,1);% p输出两个值第⼀个是a  第⼆个是b  y=ax+b
target2=polyval(squeeze(p2(:)),1:32);
clear j i wDTR2 T22
%求P值
t=1:32;
t=t';
t=[ones(32,1),t];
NCEP = NCEP';
[b1,bint1,r1,rint1,stats1]=regress(NCEP,t,0.05);
Obs = Obs';
[b2,bint2,r2,rint2,stats2]=regress(Obs,t,0.05);
%作图
clf
subplot(1,2,1)
plot(NCEP,'r','linewidth',2)
hold on
plot(Obs,'b','linewidth',2)
plot(target1,'r','linewidth',2)%趋势线
plot(target2,'b','linewidth',2)
legend({'NCEP','Observed'},'Location','Northwest')
set(gca,'xtick',[271217222732]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([-inf inf -22]);
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([133])%x轴范围锁定
box off %去掉外框
%添加斜率⽅程、R2、P值
h=text(9,-1.4,['y(NCEP)='num2str(roundn(b1(1,:),-3))'+'num2str(roundn(b1(2,:),-3))'*x''']); set(h,'FontSize',10,'FontWeight','Bold')
h=text(9,-1.6,['R^2(NCEP)='num2str(roundn(stats1(:,1),-2))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(9,-1.8,['P(NCEP)='num2str(roundn(stats1(:,3),-4))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(22,-1.4,['y(Obs)='num2str(roundn(b2(1,:),-3))'+'num2str(roundn(b2(2,:),-3))'*x''']); set(h,'FontSize',10,'FontWeight','Bold')
h=text(22,-1.6,['R^2(Obs)='num2str(roundn(stats2(:,1),-2))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(22,-1.8,['P(Obs)='num2str(roundn(stats2(:,3),-4))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(20,1.8,['Jiangxi(winter)'])%添加省份名称
set(h,'FontSize',15)
%%
%夏季
%对NCEP数据
T_mean=squeeze(nanmean(smDTR1,3));
for i=1:33
smDTR1(:,:,i)=squeeze(smDTR1(:,:,i))-T_mean;
end
T12=smDTR1;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T12(i,j,:)=nan;
end
end
end
T12=reshape(T12,[72*12833]);
T12=nanmean(T12,1);
NCEP=T12;
p1=polyfit(1:33,T12,1);% p输出两个值第⼀个是a  第⼆个是b  y=ax+b target1=polyval(squeeze(p1(:)),1:33);
clear smDTR1 T12 T_mean
%对观测数据
T_mean=squeeze(nanmean(smDTR2,3));
for i=1:33
smDTR2(:,:,i)=squeeze(smDTR2(:,:,i))-T_mean;
end
clear i j T_mean
T22=smDTR2;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T22(i,j,:)=nan;
end
end
end
T22=reshape(T22,[72*12833]);
T22=nanmean(T22,1);
Obs=T22;
p2=polyfit(1:33,T22,1);% p输出两个值第⼀个是a  第⼆个是b  y=ax+b target2=polyval(squeeze(p2(:)),1:33);
clear j i smDTR2 T22
%求P值
t=1:33;
t=t';
t=[ones(33,1),t];
NCEP = NCEP';
[b1,bint1,r1,rint1,stats1]=regress(NCEP,t,0.05);
Obs = Obs';
[b2,bint2,r2,rint2,stats2]=regress(Obs,t,0.05);
%作图
subplot(1,2,2)
plot(NCEP,'r','linewidth',2)
hold on
plot(Obs,'b','linewidth',2)
plot(target1,'r','linewidth',2)%趋势线
plot(target2,'b','linewidth',2)
legend({'NCEP','Observed'},'Location','Northwest')
set(gca,'xtick',[271217222732]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'}); axis([-inf inf -22]);
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([133])%x轴范围锁定
box off %去掉外框
%添加斜率⽅程、R2、P值
h=text(9,-1.4,['y(NCEP)='num2str(roundn(b1(1,:),-3))'+'num2str(roundn(b1(2,:),-3))'*x''']); set(h,'FontSize',10,'FontWeight','Bold')
h=text(9,-1.6,['R^2(NCEP)='num2str(roundn(stats1(:,1),-2))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(9,-1.8,['P(NCEP)='num2str(roundn(stats1(:,3),-4))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(22,-1.4,['y(Obs)='num2str(roundn(b2(1,:),-3))'+'num2str(roundn(b2(2,:),-3))'*x''']); set(h,'FontSize',10,'FontWeight','Bold')
h=text(22,-1.6,['R^2(Obs)='num2str(roundn(stats2(:,1),-2))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(22,-1.8,['P(Obs)='num2str(roundn(stats2(:,3),-4))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(18,1.8,['Jiangxi(summer)'])%添加省份名称
set(h,'FontSize',15)
得到结果如下图所⽰：
江西省冬季和夏季的DTR趋势17.2 全国平均温度的四季变化趋势
clc;clear
load('E:\数据\201507\MASK');%MASK是城市区域掩膜
load('E:\数据\201501\T1');%导⼊NCEP数据平均温度72*128*420fontweight属性bold
load('E:\数据\201501\T2');%导⼊观测数据平均温度72*128*420
%全国春季年平均温度72*128*35
for i=1:35
sT1(:,:,i)=squeeze(nanmean(T1(:,:,(3+12*(i-1)):(5+12*(i-1))),3));
sT2(:,:,i)=squeeze(nanmean(T2(:,:,(3+12*(i-1)):(5+12*(i-1))),3));
end
%%
%对NCEP数据
T_mean=squeeze(nanmean(sT1,3));
for i=1:35
sT1(:,:,i)=squeeze(sT1(:,:,i))-T_mean;
end
clear i j T_mean
T12=sT1;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T12(i,j,:)=nan;
end
end
end
T12=reshape(T12,[72*12835]);
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 j i T12
%对观测数据
T_mean=squeeze(nanmean(sT2,3));
for i=1:35
sT2(:,:,i)=squeeze(sT2(:,:,i))-T_mean;
end
clear i j T_mean
T22=sT2;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T22(i,j,:)=nan;
end
end
end
T22=reshape(T22,[72*12835]);
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 j i T22
%求P值
t=1:35;
t=t';
t=[ones(35,1),t];
NCEP = NCEP';
[b1,bint1,r1,rint1,stats1]=regress(NCEP,t,0.05);
Obs = Obs';
[b2,bint2,r2,rint2,stats2]=regress(Obs,t,0.05);
%作图
clf
subplot(1,2,1)
plot(NCEP,'r','linewidth',2)
hold on
plot(Obs,'b','linewidth',2)
plot(target1,'r','linewidth',2)%趋势线
plot(target2,'b','linewidth',2)
legend({'NCEP','Observed'},'Location','Northwest')
set(gca,'xtick',[271217222732]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([-inf inf -22]);
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([135])%x轴范围锁定为1~35
box off %去掉外框
%添加斜率⽅程、R2、P值
h=text(10,-1.5,['y(NCEP)='num2str(roundn(b1(1,:),-3))'+'num2str(roundn(b1(2,:),-3))'*x''']); set(h,'FontSize',10,'FontWeight','Bold')
h=text(10,-1.7,['R^2(NCEP)='num2str(roundn(stats1(:,1),-2))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(10,-1.9,['P(NCEP)='num2str(roundn(stats1(:,3),-4))'']);
set(h,'FontSize',10,'FontWeight','Bold')
h=text(23,-1.5,['y(Obs)='num2str(roundn(b2(1,:),-3))'+'num2str(roundn(b2(2,:),-3))'*x''']); set(h,'FontSize',10,'FontWeight','Bold')
h=text(23,-1.7,['R^2(Obs)='num2str(roundn(stats2(:,1),-2))'']);

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。