Matlab处理气象数据（五）NCEP数据的剪裁、插值和两套数据的对比

为了更直观地对比NCEP数据和观测数据在中国大陆哪些区域中上升或者下降，需要对二者进行处理。NCEP数据为全球数据，要对其进行剪裁，提取出中国区域，并将分辨率转换为与观测数据相同的0.5°×0.5°的格点数据，此处可以使用interp2函数进行插值。

首先导入NECP数据中的TemData、LonData、LatData，剪裁与插值的代码如下：

TemData=permute(TemData,[2 1 3]); %改变维数顺序，即将<192x94x420 >变为<94x192x420 >
[X,Y]=meshgrid(LonData,LatData);
[xi,yi]=meshgrid(72:0.5:135.5,18:0.5:53.5); %只要72~135.5E，18~53.5N这个范围的数据
for i=1:420
	vi(:,:,i)=interp2(X,Y,TemData(:,:,i),xi,yi); %vi为分辨率为0.5的TemData
end

此时，已经得到了两套范围和分辨率均相同的数据，但NCEP数据在中国区域之外的范围仍有数据，我们应先将其删去。

load('vi.mat');%导入NCEP数据
load('T.mat');%导入观测数据
%比较两个数组，如果数组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

用pcolor函数可以查看数据生成的图像，以1979年8月为例，分别执行

pcolor(vi(:,:,8));
shading flat;
pcolor(T(:,:,8));
shading flat;

得到处理后的两套数据范围图。（图略）

%求两套数组的年平均气温
T2=reshape(T,[72*128,420]);
Tem2=nanmean(T2,1);
Tem2=reshape(Tem2,[12,35]);
Tem2=mean(Tem2,1);

T1=reshape(vi,[72*128,420]);
Tem1=nanmean(T1,1);
Tem1=reshape(Tem1,[12,35]);
Tem1=mean(Tem1,1);

用plot函数查看两套数据的年平均气温折线图，分别执行plot(Tem1)和plot(Tem2)，得到：

两套数据的年平均气温折线图

为了更直观地展示两套数据，可以用hold on函数将两张图放在一起：

plot(Tem1,'r.-','linewidth',2); %画NCEP数据年平均气温折线图，红色实线实心点，线宽为2
hold on
plot(Tem2,'b.-','linewidth',2);
legend({'NCEP','Observed'},'Location','Northwest'); %添加图例
set(gca,'xtick',[2 7 12 17 22 27 32],'xticklabel',{'1980','1985','1990','1995','2000','2005','2010'}); %在x轴特定位置上添加标注
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

得到的图像如下所示：

两套数据年平均温度对比图

%求两套数据的差值
DT=Tem2-Tem1;
plot(DT,'b.-','linewidth',2);
hold on
plot (DT-detrend(DT),'--','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

得到下图：

两套数据年平均温度的差异

以上是平均温度的对比，最高温度和最低温度的对比与之类似。NCEP数据在NOAA网站下载最高温度（tmax.2m.mon.mean.nc）和最低温度（tmin.2m.mon.mean.nc）文件。需要注意的是，由于观测数据的最高温度和最低温度缺失2012年11、12月和2013年4、5、6、7月，所以对于最高、最低温度我们只分析1979-2011年的数据。

%提取NCEP的最高温度数据
clc,clear;
AirData  = ncread('E:\数据\空间数据\NECP\tmax.2m.mon.mean.nc','tmax'); %读入变量air
TimeData  = ncread('E:\数据\空间数据\NECP\tmax.2m.mon.mean.nc','time'); %读入变量time
LonData  = ncread('E:\数据\空间数据\NECP\tmax.2m.mon.mean.nc','lon'); %读入变量lon
LatData  = ncread('E:\数据\空间数据\NECP\tmax.2m.mon.mean.nc','lat'); %读入变量lat
AirData=squeeze(AirData);%除去size为1的维度，变为三维
TemData = AirData-273.15; %开尔文温度转换为摄氏度
TemData=permute(TemData,[2 1 3]);
[X,Y]=meshgrid(LonData,LatData);
[xi,yi]=meshgrid(72:0.5:135.5,18:0.5:53.5);%只要72~135.5E，18~53.5N这个范围的数据
for i=1:420
	vimax(:,:,i)=interp2(X,Y,TemData(:,:,i),xi,yi);%vimax为分辨率为0.5的TemData
end

————————————————————————————————

%提取NCEP的最低温度数据
clc,clear;
AirData  = ncread('E:\数据\空间数据\NECP\tmin.2m.mon.mean.nc','tmin'); %读入变量air
TimeData  = ncread('E:\数据\空间数据\NECP\tmin.2m.mon.mean.nc','time'); %读入变量time
LonData  = ncread('E:\数据\空间数据\NECP\tmin.2m.mon.mean.nc','lon'); %读入变量lon
LatData  = ncread('E:\数据\空间数据\NECP\tmin.2m.mon.mean.nc','lat'); %读入变量lat
AirData=squeeze(AirData);%除去size为1的维度，变为三维
TemData = AirData-273.15; %开尔文温度转换为摄氏度
TemData=permute(TemData,[2 1 3]);
[X,Y]=meshgrid(LonData,LatData);
[xi,yi]=meshgrid(72:0.5:135.5,18:0.5:53.5);%只要72~135.5E，18~53.5N这个范围的数据
for i=1:420
	vimin(:,:,i)=interp2(X,Y,TemData(:,:,i),xi,yi);%vimin为分辨率为0.5的TemData
end

---

%提取观测数据最高温度数据
clc,clear
for iyear=1979:2011
	for imonth=1:12
    	if(imonth>0&&imonth<=9)
		Tpath=strcat('E:\数据\MAX\SURF_CLI_CHN_TEM_MON_GRID_0.5-MAX-',num2str(iyear),'0',num2str(imonth),'.txt'); 
   		 else
		Tpath=strcat('E:\数据\MAX\SURF_CLI_CHN_TEM_MON_GRID_0.5-MAX-',num2str(iyear),'',num2str(imonth),'.txt');     
  		end
		Tmax(imonth,iyear-1978,:,:)=textread(Tpath);
	end
end

for i=1:12
	for j=1:33
		for k=1:72
			for l=1:128
				if Tmax(i,j,k,l)==-99.0;
				Tmax(i,j,k,l)=nan;
				end
			end
		end
	end
end

Tmax=reshape(Tmax,[396,72*128]);
Tem=nanmean(Tmax,2);
Tem=reshape(Tem,[12,33]);
Tem2=mean(Tem,1);
%plot(Tem2);
Tmax=reshape(Tmax,[396,72,128]);
Tmax=permute(Tmax,[2 3 1]); %使之成为72*128*396型的矩阵
%使图像颠倒，数据从最左下角向右上角读起
for i=1:396
	Tmax (:,:,i)=flipud(squeeze(Tmax (:,:,i)));
end

---

%提取观测数据最低温度数据
clc,clear
for iyear=1979:2011
	for imonth=1:12
   		if(imonth>0&&imonth<=9)
		Tpath=strcat('E:\数据\MIN\SURF_CLI_CHN_TEM_MON_GRID_0.5-MIN-',num2str(iyear),'0',num2str(imonth),'.txt'); 
		else
		Tpath=strcat('E:\数据\MIN\SURF_CLI_CHN_TEM_MON_GRID_0.5-MIN-',num2str(iyear),'',num2str(imonth),'.txt');     
		end
		Tmin(imonth,iyear-1978,:,:)=textread(Tpath);
	end
end
 
for i=1:12
	for j=1:33
		for k=1:72
			for l=1:128
				if Tmin(i,j,k,l)==-99.0;
				Tmin(i,j,k,l)=nan;
				end
			end
		end
	end
end
 
Tmin=reshape(Tmin,[396,72*128]);
Tem=nanmean(Tmin,2);
Tem=reshape(Tem,[12,33]);
Tem2=mean(Tem,1);
%plot(Tem2);
Tmin=reshape(Tmin,[396,72,128]);
Tmin=permute(Tmin,[2 3 1]); %使之成为72*128*396型的矩阵
%使图像颠倒，数据从最左下角向右上角读起
for i=1:396
	Tmin(:,:,i)=flipud(squeeze(Tmin(:,:,i)));
end

至此，处理后的NCEP数据和观测数据平均温度、最高温度和最低温度分别保存为vi.mat、vimax.mat、vimin.mat，T.mat、Tmax.mat、Tmin.mat。

---

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