欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识！
您现在的位置是： 首页

Matlab处理气象数据（十六）城市与非城市区域的对比

程序员文章站 2022-07-14 11:54:35
...

利用2010城市边界数据，建立掩膜。从全国数据中抠出城市区域的数据，然后反向抠出非城市区域的数据，最后对比两个区域趋势变化的不同。
Matlab处理气象数据（十六）城市与非城市区域的对比

栅格化的城市区域掩膜

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

得到结果如下图所示：
Matlab处理气象数据（十六）城市与非城市区域的对比

城市区域的异常和趋势线

Matlab处理气象数据（十六）城市与非城市区域的对比

非城市区域的异常和趋势线

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

结果如下图所示：
Matlab处理气象数据（十六）城市与非城市区域的对比

城市区域DTR的异常和趋势

Matlab处理气象数据（十六）城市与非城市区域的对比

非城市区域DTR的异常和趋势

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处理气象数据（十六）城市与非城市区域的对比

城市和非城市区域DTR异常和趋势

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

相关标签： Matlab 气象 matlab

上一篇： K-近邻算法（一）

下一篇： 机器学习--k近邻算法

推荐阅读