Matlab处理气象数据（六）数据的面积加权

由于地理投影导致导致每个像元实际地面面积不同，越靠近北极实际面积越小，越靠近赤道实际面积越大，如果不进行面积加权就简单平均，会导致温度较实际温度偏低。用以下公式计算地理投影每个像元的面积：double pi = 3.1415926;double R = 6371007.181; //units - mdouble area = (pi/180.0)*R*R*fabs(sin(dLat1/180.0

岩鸣杨子

10444人浏览 · 2020-05-10 18:23:25

岩鸣杨子 · 2020-05-10 18:23:25 发布

由于地理投影导致导致每个像元实际地面面积不同，越靠近北极实际面积越小，越靠近赤道实际面积越大，如果不进行面积加权就简单平均，会导致温度较实际温度偏低。

用以下公式计算地理投影每个像元的面积：

double pi = 3.1415926;
double R = 6371007.181; //units - m						
double area = (pi/180.0)*R*R*fabs(sin(dLat1/180.0*pi) - sin(dLat2/180.0*pi)) * fabs(dLon1 - dLon2);

其中，dLat1, dLat2是像元上下边界（北南）的纬度，dLon1, dLon2是像元左右边界的经度。

下面分别对观测数据和NCEP数据的平均、最高、最低温度进行面积加权平均：

%对观测数据平均温度进行面积加权
load('T.mat'); %导入观测数据72*128*420

%计算实际面积
R =  6371007.181;
dLat1=18.5:0.5:54;
dLat2=18:0.5:53.5;
dLon1=72:0.5:135.5;
dLon2=72.5:0.5:136;

for i=1:72
	for j=1:128
		area(i,j) = (pi/180.0)*R*R*abs(sin(dLat1(i)/180.0*pi) - sin(dLat2(i)/180.0*pi)) * abs(dLon1(j) - dLon2(j));
	end
end

%将T为NAN时area对应位置上也设为NAN
for i=1:72
	for j=1:128
		if isnan(T(i,j,1))==1
			area(i,j)=nan;
		end
	end
end
a = nansum(area(:));%求像元的总面积

%求面积的权值W
for i=1:72
	for j=1:128
		W(i,j)= area(i,j)/a;
	end
end

%计算加权后的温度
WT=zeros(1,420);
for t=1:420
	WT(t)=0;
	for i=1:72
		for j=1:128
			if isnan(T(i,j,t))==1;
				continue
			else
				WT(t)=WT(t)+T(i,j,t).*W(i,j);
			end
		end
	end
end
WT2=reshape(WT,[12,35]);
Tem2=mean(WT2,1);
plot(Tem2);

%对观测数据最高温度进行面积加权
load('Tmax.mat'); %导入观测数据72*128*396

%计算实际面积
R =  6371007.181;
dLat1=18.5:0.5:54;
dLat2=18:0.5:53.5;
dLon1=72:0.5:135.5;
dLon2=72.5:0.5:136;

for i=1:72
	for j=1:128
		area(i,j) = (pi/180.0)*R*R*abs(sin(dLat1(i)/180.0*pi) - sin(dLat2(i)/180.0*pi)) * abs(dLon1(j) - dLon2(j));
	end
end

%将Tmax为NAN时area对应位置上也设为NAN
for i=1:72
	for j=1:128
		if isnan(Tmax(i,j,1))==1
			area(i,j)=nan;
		end
	end
end
a = nansum(area(:));%求像元的总面积

%求面积的权值W
for i=1:72
	for j=1:128
		W(i,j)= area(i,j)/a;
	end
end

%计算加权后的温度
WT=zeros(1,396);
for t=1:396
	WT(t)=0;
	for i=1:72
		for j=1:128
			if isnan(Tmax(i,j,t))==1;
				continue
			else
				WT(t)=WT(t)+Tmax(i,j,t).*W(i,j);
			end
		end
	end
end
WTmax2=reshape(WT,[12,33]);
Temmax2=mean(WTmax2,1);
plot(Temmax2);

%对观测数据最低温度进行面积加权
load('Tmin.mat'); %导入观测数据72*128*396

%计算实际面积
R =  6371007.181;
dLat1=18.5:0.5:54;
dLat2=18:0.5:53.5;
dLon1=72:0.5:135.5;
dLon2=72.5:0.5:136;

for i=1:72
	for j=1:128
		area(i,j) = (pi/180.0)*R*R*abs(sin(dLat1(i)/180.0*pi) - sin(dLat2(i)/180.0*pi)) * abs(dLon1(j) - dLon2(j));
	end
end

%将Tmin为NAN时area对应位置上也设为NAN
for i=1:72
	for j=1:128
		if isnan(Tmin(i,j,1))==1
			area(i,j)=nan;
		end
	end
end
a = nansum(area(:));%求像元的总面积

%求面积的权值W
for i=1:72
	for j=1:128
		W(i,j)= area(i,j)/a;
	end
end

%计算加权后的温度
WT=zeros(1,396);
	for t=1:396
		WT(t)=0;
		for i=1:72
			for j=1:128
				if isnan(Tmin(i,j,t))==1;
					continue
				else
					WT(t)=WT(t)+Tmin(i,j,t).*W(i,j);
			end
		end
	end
end
WTmin2=reshape(WT,[12,33]);
Temmin2=mean(WTmin2,1);
plot(Temmin 2);

%对NCEP数据平均温度进行面积加权
load('vi.mat');%导入NCEP数据72*128*420
load('area.mat');%导入面积数据72*128
load('T.mat');%导入观测数据72*128*420

%比较两个数组，如果数组T为NAN时，数组vi相对应位置也设为NAN。
for i=1:420
	for j=1:128
		for k=1:72
			if isnan(T(k,j,i))==1
				vi(k,j,i)=nan;
			end
		end
	end
end

%将vi为NAN时area对应位置上也设为NAN
for i=1:72
	for j=1:128
		if isnan(vi(i,j,1))==1
			area(i,j)=nan;
		end
	end
end
a = nansum(area(:)); %求像元的总面积

%求面积的权值W
for i=1:72
	for j=1:128
		W(i,j)= area(i,j)/a;
	end
end

%计算加权后的温度
WT1=zeros(1,420);
for t=1:420
	WT1(t)=0;
	for i=1:72
		for j=1:128
			if isnan(vi(i,j,t))==1;
				continue
			else
				WT1(t)=WT1(t)+vi(i,j,t).*W(i,j);
			end
		end
	end
end
Tem1=reshape(WT1,[12,35]);
Tem1=mean(Tem1,1);
plot(Tem1);

%对NCEP数据最高温度进行面积加权
load('vimax.mat');%导入NCEP数据72*128*420
load('area.mat');%导入面积数据72*128
load('Tmax.mat');%导入观测数据72*128*420

%比较两个数组，如果数组Tmax为NAN时，数组vimax相对应位置也设为NAN。
for i=1:396
	for j=1:128
		for k=1:72
			if isnan(Tmax(k,j,i))==1
				vimax(k,j,i)=nan;
			end
		end
	end
end

%将vimax为NAN时area对应位置上也设为NAN
for i=1:72
	for j=1:128
		if isnan(vimax(i,j,1))==1
			area(i,j)=nan;
		end
	end
end
a = nansum(area(:));%求像元的总面积

%求面积的权值W
for i=1:72
	for j=1:128
		W(i,j)= area(i,j)/a;
	end
end

%计算加权后的温度
WTmax1=zeros(1,420);
for t=1:420
	WTmax1(t)=0;
	for i=1:72
		for j=1:128
			if isnan(vimax(i,j,t))==1;
				continue
			else
				WTmax1(t)=WTmax1(t)+vimax(i,j,t).*W(i,j);
			end
		end
	end
end
Temmax1=reshape(WTmax1,[12,35]);
Temmax1=mean(Temmax1,1);
plot(Temmax1);

%对NCEP数据最低温度进行面积加权
load('vimin.mat');%导入NCEP数据72*128*420
load('area.mat');%导入面积数据72*128
load('Tmin.mat');%导入观测数据72*128*420

%比较两个数组，如果数组Tmin为NAN时，数组vimin相对应位置也设为NAN。
for i=1:396
	for j=1:128
		for k=1:72
			if isnan(Tmin(k,j,i))==1
				vimin(k,j,i)=nan;
			end
		end
	end
end

%将vimin为NAN时area对应位置上也设为NAN
for i=1:72
	for j=1:128
		if isnan(vimin(i,j,1))==1
			area(i,j)=nan;
		end
	end
end
a = nansum(area(:));%求像元的总面积

%求面积的权值W
for i=1:72
	for j=1:128
		W(i,j)= area(i,j)/a;
	end
end

%计算加权后的温度
WTmin1=zeros(1,420);
for t=1:420
	WTmin1(t)=0;
	for i=1:72
		for j=1:128
			if isnan(vimin(i,j,t))==1;
				continue
			else
				WTmin1(t)=WTmin1(t)+vimin(i,j,t).*W(i,j);
			end
		end
	end
end
Temmin1=reshape(WTmin1,[12,35]);
Temmin1=mean(Temmin1,1);
plot(Temmin1);

至此，进行面积加权处理后的NCEP数据和观测数据平均温度、最高温度和最低温度分别保存为Tem1 .mat、Temmax1.mat、Temmin1.mat，Tem2.mat、Temmax2.mat、Temmin2.mat。

%面积加权后的两套数据平均、最高、最低温度对比
DT = Tem2 - Tem1;
DTmax = Temmax2 - Temmax1;
DTmin = Temmin2 - Temmin1;
plot(DT,'r.-','linewidth',2);
hold on
plot(DTmax,'b.-','linewidth',2);
plot(DTmin,'g.-','linewidth',2);
legend({'DT','DTmax','DTmin'},'Location','Northwest');%添加图例
plot (DTmax-detrend(DTmax),'b--','linewidth',2);
plot (DTmin-detrend(DTmin),'g--','linewidth',2);
plot (DT-detrend(DT),'r--','linewidth',2);
set(gca,'xtick',[2 7 12 17 22 27 32],'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'})
set(gca, 'FontSize',10,'FontWeight','Bold','tickdir','out') %设置标注为10号字、加粗、标记线向外
h=xlabel('Year'); %设置x轴名称
set(h, 'FontSize',10,'FontWeight','Bold')
h=ylabel('Temperarure(\circC)'); %设置y轴名称
set(h, 'FontSize',10,'FontWeight','Bold') 
xlim([1 35])%x轴范围锁定为1~35
box off %去掉外框
hold off