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*128 32]);
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 j=1:128
if isnan(MASK(i,j))==1
T22(i,j,:)=nan;
end
end
end
T22=reshape(T22,[72*128 32]);
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',[2 7 12 17 22 27 32]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([ -inf inf -2 2]);
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 33])%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*128 33]);
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*128 33]);
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',[2 7 12 17 22 27 32]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([ -inf inf -2 2]);
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 33])%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)
得到结果如下图所示:
17.2 全国平均温度的四季变化趋势
clc;clear
load('E:\数据\201507\MASK');%MASK是城市区域掩膜
load('E:\数据\201501\T1');%导入NCEP数据平均温度 72*128*420
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*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 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*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 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',[2 7 12 17 22 27 32]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([ -inf inf -2 2]);
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 %去掉外框
%添加斜率方程、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)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(23,-1.9,['P(Obs)=' num2str(roundn(stats2(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(26,1.8,['spring(a)']) %添加名称
set(h, 'FontSize',15)
%可以用 delete(h) 来撤销句柄
%%
%对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
%抠非城市的
for i=1:72
for j=1:128
if isnan(sT1(i,j,:))==1
sT1(i,j,:)=0;
end
end
end
for i=1:72
for j=1:128
if isnan(T12(i,j,:))==1
T12(i,j,:)=0;
end
end
end
T13=sT1-T12;
for i=1:72
for j=1:128
if T13(i,j,:)==0
T13(i,j,:)=nan;
end
end
end
T13=reshape(T13,[72*128 35]);
T13=nanmean(T13,1);
NCEP=T13;
p1=polyfit(1:35,T13,1); % p输出两个值第一个是a 第二个是b y=ax+b
target1=polyval(squeeze(p1(:)),1:35);
clear j i sT1 T13
%对观测数据
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
%抠非城市的
for i=1:72
for j=1:128
if isnan(sT2(i,j,:))==1
sT2(i,j,:)=0;
end
end
end
for i=1:72
for j=1:128
if isnan(T22(i,j,:))==1
T22(i,j,:)=0;
end
end
end
T23=sT2-T22;
for i=1:72
for j=1:128
if T23(i,j,:)==0
T23(i,j,:)=nan;
end
end
end
T23=reshape(T23,[72*128 35]);
T23=nanmean(T23,1);
Obs=T23;
p2=polyfit(1:35,T23,1); % p输出两个值第一个是a 第二个是b y=ax+b
target2=polyval(squeeze(p2(:)),1:35);
clear j i sT2 T23
%求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);
%作图
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',[2 7 12 17 22 27 32]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([ -inf inf -2 2]);
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 %去掉外框
%添加斜率方程、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)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(23,-1.9,['P(Obs)=' num2str(roundn(stats2(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(26,1.8,['spring(b)']) %添加名称
set(h, 'FontSize',15)
%可以用 delete(h) 来撤销句柄
类似地,可以得到全国城市区域和非城市区域最高温度和最低温度的四季变化趋势,此处不再赘述。
17.3 全国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*33
for i=1:33
sDTR1(:,:,i)=squeeze(nanmean(DTR1(:,:,(3+12*(i-1)):(5+12*(i-1))),3));
sDTR2(:,:,i)=squeeze(nanmean(DTR2(:,:,(3+12*(i-1)):(5+12*(i-1))),3));
end
%%
%对NCEP数据
T_mean=squeeze(nanmean(sDTR1,3));
for i=1:33
sDTR1(:,:,i)=squeeze(sDTR1(:,:,i))-T_mean;
end
clear i j T_mean
T12=sDTR1;%删除边界区域外的值
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*128 33]);
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 j i T12
%对观测数据
T_mean=squeeze(nanmean(sDTR2,3));
for i=1:33
sDTR2(:,:,i)=squeeze(sDTR2(:,:,i))-T_mean;
end
clear i j T_mean
T22=sDTR2;%删除边界区域外的值
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*128 33]);
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 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);
%作图
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',[2 7 12 17 22 27 32]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([ -inf inf -2 2]);
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 33])%x轴范围锁定为1~33
box off %去掉外框
%添加斜率方程、R2、P值
h=text(8,-1.5,['y(NCEP)=' num2str(roundn(b1(1,:),-3)) '+' num2str(roundn(b1(2,:),-3)) '*x' '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(8,-1.7,['R^2(NCEP)=' num2str(roundn(stats1(:,1),-2)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(8,-1.9,['P(NCEP)=' num2str(roundn(stats1(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(21,-1.5,['y(Obs)=' num2str(roundn(b2(1,:),-3)) '+' num2str(roundn(b2(2,:),-3)) '*x' '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(21,-1.7,['R^2(Obs)=' num2str(roundn(stats2(:,1),-2)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(21,-1.9,['P(Obs)=' num2str(roundn(stats2(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(26,1.8,['spring(a)']) %添加名称
set(h, 'FontSize',15)
%可以用 delete(h) 来撤销句柄
%%
%对NCEP数据
T_mean=squeeze(nanmean(sDTR1,3));
for i=1:33
sDTR1(:,:,i)=squeeze(sDTR1(:,:,i))-T_mean;
end
clear i j T_mean
T12=sDTR1;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T12(i,j,:)=nan;
end
end
end
%抠非城市的
for i=1:72
for j=1:128
if isnan(sDTR1(i,j,:))==1
sDTR1(i,j,:)=0;
end
end
end
for i=1:72
for j=1:128
if isnan(T12(i,j,:))==1
T12(i,j,:)=0;
end
end
end
T13=sDTR1-T12;
for i=1:72
for j=1:128
if T13(i,j,:)==0
T13(i,j,:)=nan;
end
end
end
T13=reshape(T13,[72*128 33]);
T13=nanmean(T13,1);
NCEP=T13;
p1=polyfit(1:33,T13,1); % p输出两个值第一个是a 第二个是b y=ax+b
target1=polyval(squeeze(p1(:)),1:33);
clear j i sDTR1 T13
%对观测数据
T_mean=squeeze(nanmean(sDTR2,3));
for i=1:33
sDTR2(:,:,i)=squeeze(sDTR2(:,:,i))-T_mean;
end
clear i j T_mean
T22=sDTR2;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T22(i,j,:)=nan;
end
end
end
%抠非城市的
for i=1:72
for j=1:128
if isnan(sDTR2(i,j,:))==1
sDTR2(i,j,:)=0;
end
end
end
for i=1:72
for j=1:128
if isnan(T22(i,j,:))==1
T22(i,j,:)=0;
end
end
end
T23=sDTR2-T22;
for i=1:72
for j=1:128
if T23(i,j,:)==0
T23(i,j,:)=nan;
end
end
end
T23=reshape(T23,[72*128 33]);
T23=nanmean(T23,1);
Obs=T23;
p2=polyfit(1:33,T23,1); % p输出两个值第一个是a 第二个是b y=ax+b
target2=polyval(squeeze(p2(:)),1:33);
clear j i sDTR2 T23
%求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',[2 7 12 17 22 27 32]);
set(gca,'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'});
axis([ -inf inf -2 2]);
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 33])%x轴范围锁定为1~33
box off %去掉外框
%添加斜率方程、R2、P值
h=text(8,-1.5,['y(NCEP)=' num2str(roundn(b1(1,:),-3)) '+' num2str(roundn(b1(2,:),-3)) '*x' '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(8,-1.7,['R^2(NCEP)=' num2str(roundn(stats1(:,1),-2)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(8,-1.9,['P(NCEP)=' num2str(roundn(stats1(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(21,-1.5,['y(Obs)=' num2str(roundn(b2(1,:),-3)) '+' num2str(roundn(b2(2,:),-3)) '*x' '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(21,-1.7,['R^2(Obs)=' num2str(roundn(stats2(:,1),-2)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(21,-1.9,['P(Obs)=' num2str(roundn(stats2(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(26,1.8,['spring(b)']) %添加名称
set(h, 'FontSize',15)
%可以用 delete(h) 来撤销句柄
相关链接:
Matlab处理气象数据——目录