利用2010城市边界数据，建立掩膜。从全国数据中抠出城市区域的数据，然后反向抠出非城市区域的数据，最后对比两个区域趋势变化的不同。
16.1 抠出城市和非城市的数据
clc;clear
%扣NCEP的城市和非城市区域数据
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:\数据\边界\2010urban_geographic\2010urban_geographic.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;%这里保存MASK数据，留到后面用
T12=T1;%删除边界区域外的值
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
T12(i,j,:)=nan;%T12即城市区域的数据
end
end
end
%抠非城市的
for i=1:72
for j=1:128
if isnan(T1(i,j,:))==1
T1(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=T1-T12;
for i=1:72
for j=1:128
if T13(i,j,:)==0
T13(i,j,:)=nan;%T13即非城市区域的数据
end
end
end
%扣观测数据的城市和非城市区域数据
lon=72:0.5:135.5;
lat=18:0.5:53.5;
[x y]=meshgrid(lon,lat);
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:\数据\边界\2010urban_geographic\2010urban_geographic.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
%抠非城市的
for i=1:72
for j=1:128
if isnan(T2(i,j,:))==1
T2(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=T2-T22;
for i=1:72
for j=1:128
if T23(i,j,:)==0
T23(i,j,:)=nan;
end
end
end
16.2 城市和非城市数据的异常和趋势
%城市区域的距平和趋势线
T12=reshape(T12,[72*128,35]);
T12=nanmean(T12,1);
T22=reshape(T22,[72*128,35]);
T22=nanmean(T22,1);
m1=mean(T12); %求Tem1的平均值
m2=mean(T22);
a1=T12-m1; %求Tem1的距平
a2=T22-m2;
plot(a1,'r.-','linewidth',2);%画NCEP数据年平均气温折线图，红色实线实心点
hold on
plot(a2,'b.-','linewidth',2);
plot(a1-detrend(a1),'r-','linewidth',2);
plot(a2-detrend(a2),'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 -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 %去掉外框
hold off
%非城市区域的距平和趋势线
T13=reshape(T13,[72*128,35]);
T13=nanmean(T13,1);
T23=reshape(T23,[72*128,35]);
T23=nanmean(T23,1);
m1=mean(T13); %求Tem1的平均值
m2=mean(T23);
a1=T13-m1; %求Tem1的距平
a2=T23-m2;
plot(a1,'r.-','linewidth',2);%画NCEP数据年平均气温折线图，红色实线实心点
hold on
plot(a2,'b.-','linewidth',2);
plot(a1-detrend(a1),'r-','linewidth',2);
plot(a2-detrend(a2),'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 -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 %去掉外框
hold off
得到结果如下图所示：
16.3 计算城市和非城市的数据的DTR
%计算NCEP数据城市和非城市区域的DTR
load('E:\数据\201507\MASK');
load('E:\数据\201501\Tmax1');
load('E:\数据\201501\Tmin1');
DTR1=Tmax1-Tmin1;
T10=reshape(DTR1,[72 128 12 33]);
T11=reshape(DTR1,[72 128 12*33]);
T10=squeeze(nanmean(T10,3));
T_mean=squeeze(nanmean(T11,3));
clear DTR1
for i=1:33
DTR1(:,:,i)=squeeze(T10(:,:,i))-T_mean;
end
clear i j T_mean T11 T10 Tmax1 Tmin1
DTR12=DTR1;
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
DTR12(i,j,:)=nan;
end
end
end
%save DTR12;%NCEP数据城市区域的DTR
for i=1:72
for j=1:128
if isnan(DTR1(i,j,:))==1
DTR1(i,j,:)=0;
end
end
end
for i=1:72
for j=1:128
if isnan(DTR12(i,j,:))==1
DTR12(i,j,:)=0;
end
end
end
DTR13=DTR1-DTR12;
for i=1:72
for j=1:128
if DTR13(i,j,:)==0
DTR13(i,j,:)=nan;
end
end
end
%save DTR13;%NCEP数据非城市区域的DTR
%计算观测数据数据城市和非城市区域的DTR
load('E:\数据\201507\MASK');
load('E:\数据\201501\Tmax2');
load('E:\数据\201501\Tmin2');
DTR2=Tmax2-Tmin2;
T10=reshape(DTR2,[72 128 12 33]);
T11=reshape(DTR2,[72 128 12*33]);
T10=squeeze(nanmean(T10,3));
T_mean=squeeze(nanmean(T11,3));
clear DTR2
for i=1:33
DTR2(:,:,i)=squeeze(T10(:,:,i))-T_mean;
end
clear i j T_mean T11 T10 Tmax2 Tmin2
DTR22=DTR2;
for i=1:72
for j=1:128
if isnan(MASK(i,j))==1
DTR22(i,j,:)=nan;
end
end
end
%save DTR22;%观测数据城市区域的DTR
for i=1:72
for j=1:128
if isnan(DTR2(i,j,:))==1
DTR2(i,j,:)=0;
end
end
end
for i=1:72
for j=1:128
if isnan(DTR22(i,j,:))==1
DTR22(i,j,:)=0;
end
end
end
DTR23=DTR2-DTR22;
for i=1:72
for j=1:128
if DTR23(i,j,:)==0
DTR23(i,j,:)=nan;
end
end
end
%save DTR23;%观测数据非城市区域的DTR
16.4 城市和非城市的数据的DTR趋势的异常和趋势线
%城市区域
load('DTR12.mat')
load('DTR22.mat')
DTR12=reshape(DTR12,[72*128,33]);
DTR12=nanmean(DTR12,1);
DTR22=reshape(DTR22,[72*128,33]);
DTR22=nanmean(DTR22,1);
m1=mean(DTR12);
m2=mean(DTR22);
a1=DTR12-m1;
a2=DTR22-m2;
plot(a1,'r.-','linewidth',2);%画NCEP数据DTR折线图，红色实线实心点
hold on
plot(a2,'b.-','linewidth',2);
plot(a1-detrend(a1),'r-','linewidth',2);
plot(a2-detrend(a2),'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 -1 1]);
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 %去掉外框
hold off
%非城市区域
load('DTR13.mat')
load('DTR23.mat')
DTR13=reshape(DTR13,[72*128,33]);
DTR13=nanmean(DTR13,1);
DTR23=reshape(DTR23,[72*128,33]);
DTR23=nanmean(DTR23,1);
m1=mean(DTR13);
m2=mean(DTR23);
a1=DTR13-m1;
a2=DTR23-m2;
plot(a1,'r.-','linewidth',2);%画NCEP数据DTR折线图，红色实线实心点
hold on
plot(a2,'b.-','linewidth',2);
plot(a1-detrend(a1),'r-','linewidth',2);
plot(a2-detrend(a2),'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 -1 1]);
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 %去掉外框
hold off
结果如下图所示：
16.5 城市和非城市的DTR的R2和P值
在上面结果的基础上，增加P值和R2值的计算，并添加斜率方程。
clc;clear
%%抠NCEP数据
%Tmax1、Tmax2、Tmin1、Tmin2：抠出来的两套中国数据的最高最低值，格式为72*128*396
%MASK是城市区域掩膜
load('E:\数据\201507\MASK');
load('E:\数据\201501\Tmax1');
load('E:\数据\201501\Tmin1');
DTR1=Tmax1-Tmin1;
T10=reshape(DTR1,[72 128 12 33]);
T11=reshape(DTR1,[72 128 12*33]);
T10=squeeze(nanmean(T10,3));
T_mean=squeeze(nanmean(T11,3));
clear DTR1
for i=1:33
DTR1(:,:,i)=squeeze(T10(:,:,i))-T_mean;
end
clear i j T_mean T11 T10 Tmax1 Tmin1
T12=DTR1;%删除边界区域外的值
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 DTR1 T12
%%
%抠观测数据
%Tmax1、Tmax2、Tmin1、Tmin2：抠出来的两套中国数据的最高最低值，格式为72*128*396
load('E:\数据\201501\Tmax2');
load('E:\数据\201501\Tmin2');
DTR2=Tmax2-Tmin2;
T10=reshape(DTR2,[72 128 12 33]);
T11=reshape(DTR2,[72 128 12*33]);
T10=squeeze(nanmean(T10,3));
T_mean=squeeze(nanmean(T11,3));
clear DTR2
for i=1:33
DTR2(:,:,i)=squeeze(T10(:,:,i))-T_mean;
end
clear i j T_mean T11 T10 Tmax2 Tmin2
T22=DTR2;%删除边界区域外的值
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 DTR2 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);
%%
%作图
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 -1 1]);
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~35
box off %去掉外框
%添加斜率方程、R2、P值
h=text(10,-0.7,['y(NCEP)=' num2str(roundn(b1(1,:),-3)) '+' num2str(roundn(b1(2,:),-3)) '*x' '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(10,-0.8,['R^2(NCEP)=' num2str(roundn(stats1(:,1),-2)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(10,-0.9,['P(NCEP)=' num2str(roundn(stats1(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(23,-0.7,['y(Obs)=' num2str(roundn(b2(1,:),-3)) '+' num2str(roundn(b2(2,:),-3)) '*x' '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(23,-0.8,['R^2(Obs)=' num2str(roundn(stats2(:,1),-2)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(23,-0.9,['P(Obs)=' num2str(roundn(stats2(:,3),-4)) '']);
set(h, 'FontSize',10,'FontWeight','Bold')
h=text(31,0.9,['(a)']) %添加省份名称
set(h, 'FontSize',15)
%以上代码可以用 delete(h) 来撤销句柄
相关链接：
Matlab处理气象数据——目录
563
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
