观测数据为中国气象数据网上获得的中国地面温度月值0.5°×0.5°格点数据集,时间范围是1961年1月至2013年12月。这套数据为txt格式,包含头文件。头文件信息为:
NCOLS 128
NROWS 72
XLLCORNER 72.000000000000
YLLCORNER 18.000000000000
CELLSIZE 0.50000000000000
NODATA_VALUE -99.00000
从中可以得知,有128列、72行,横轴坐标从72度开始,纵轴坐标从18开始,每个格点为0.5度,缺失值默认填充为-99.0。由此可以推知,本数据坐标范围是[72 136]*[18 54]。
由于数据带有头文件,不能直接导入matlab中进行计算,所以应先去除头文件,再将缺失值全部替换为NAN。批量替换代码如下:
clc,clear
for iyear=1979:2013
for imonth=1:12
if(imonth>0&&imonth<=9)
Tpath=strcat('E:\数据\MEAN\SURF_CLI_CHN_TEM_MON_GRID_0.5-MEAN-',num2str(iyear),'0',num2str(imonth),'.txt');
else
Tpath=strcat('E:\数据\MEAN\SURF_CLI_CHN_TEM_MON_GRID_0.5-MEAN-',num2str(iyear),'',num2str(imonth),'.txt');
end
T(imonth,iyear-1978,:,:)=textread(Tpath);
end
end
for i=1:12
for j=1:35
for k=1:72
for l=1:128
if T(i,j,k,l)==-99.0;
T(i,j,k,l)=nan; %将缺失值替换为NAN
end
end
end
end
end
此时,T为4-D数组,其格式为<month,year,lon,lat>
。
T=reshape(T,[420,72,128]);
T=permute(T,[2 3 1]);%使之成为72*128*420型的矩阵
注意:这套观测数据是从左上角开始读取的,如果不作调整直接出图,则图像会是颠倒的。这里使用flipud
函数,对矩阵实现上下翻转。
%使图像颠倒,数据从最左下角向右上角读起
for i=1:420
T(:,:,i)=flipud(squeeze(T(:,:,i)));
end
要计算中国区域的年平均值,则可先求出每个月的平均值,再对每12个月求平均。
T=reshape(T,[420,72*128]); %<420x9216 double>
Tem=nanmean(T,2); %420个月的平均值<420x1 double>
Tem=reshape(Tem,[12,35]);%<12x35 double>
Tem2=mean(Tem,1); % %将观测数据年均值命名为Tem2<1x35 double>
plot(Tem2);
得到下图:
对比之前用NCEP数据得到的年平均气温变化图,可以看出两者都呈现上升的趋势,但观测数据比NCEP数据温度更高。
注:此处还未对数据进行加权处理。
相关链接:
Matlab处理气象数据——目录