matlab风速数据book,基于Matlab的中国地面气候资料日值数据集(V3.0)快速处理

本次主要是对中国气象数据网的V3.0地面气象观测数据集进行快速处理，将原先分散的txt各个要素整合到一个列表中，一行代表一个站点。直接从全国的数据集中来提取出来自己想要的站点数据集，该过程也包含了一年缺测数据低于30天的站点数据线性内插，同时参考数据说明对异常值进行了处理。本次以新疆省内的气象站点为例子，从全国V3.0数据集中提取新疆站点。

% author email: yinlichang3064@163.com

station=xlsread('D:\xj\新疆.xlsx');% 一列新疆的站点数据

%

for year=2000:2017

datasum=[];

for mon=1:12

if mon<10

filename_pre=['F:\Day1951-2012yuanshi\DAY_PRE195101201612\SURF_CLI_CHN_MUL_DAY-PRE-13011-',int2str(year),'0',int2str(mon),'.TXT'];

filename_temp=['F:\Day1951-2012yuanshi\DAY_TEM195101201612\SURF_CLI_CHN_MUL_DAY-TEM-12001-',int2str(year),'0',int2str(mon),'.TXT'];

filename_rhu=['F:\Day1951-2012yuanshi\DAY_RHU195101201612\SURF_CLI_CHN_MUL_DAY-RHU-13003-',int2str(year),'0',int2str(mon),'.TXT'];

filename_SSD=['F:\Day1951-2012yuanshi\DAY_SSD195101201612\SURF_CLI_CHN_MUL_DAY-SSD-14032-',int2str(year),'0',int2str(mon),'.TXT'];

filename_WIN=['F:\Day1951-2012yuanshi\DAY_WIN195101201612\SURF_CLI_CHN_MUL_DAY-WIN-11002-',int2str(year),'0',int2str(mon),'.TXT'];

filename_prs=['F:\Day1951-2012yuanshi\DAY_PRS195101201612\SURF_CLI_CHN_MUL_DAY-PRS-10004-',int2str(year),'0',int2str(mon),'.TXT'];

else

filename_pre=['F:\Day1951-2012yuanshi\DAY_PRE195101201612\SURF_CLI_CHN_MUL_DAY-PRE-13011-',int2str(year),int2str(mon),'.TXT'];

filename_temp=['F:\Day1951-2012yuanshi\DAY_TEM195101201612\SURF_CLI_CHN_MUL_DAY-TEM-12001-',int2str(year),int2str(mon),'.TXT'];

filename_rhu=['F:\Day1951-2012yuanshi\DAY_RHU195101201612\SURF_CLI_CHN_MUL_DAY-RHU-13003-',int2str(year),int2str(mon),'.TXT'];

filename_SSD=['F:\Day1951-2012yuanshi\DAY_SSD195101201612\SURF_CLI_CHN_MUL_DAY-SSD-14032-',int2str(year),int2str(mon),'.TXT'];

filename_WIN=['F:\Day1951-2012yuanshi\DAY_WIN195101201612\SURF_CLI_CHN_MUL_DAY-WIN-11002-',int2str(year),int2str(mon),'.TXT'];

filename_prs=['F:\Day1951-2012yuanshi\DAY_PRS195101201612\SURF_CLI_CHN_MUL_DAY-PRS-10004-',int2str(year),int2str(mon),'.TXT'];

end

data_prec=importdata(filename_pre);

data_temp=importdata(filename_temp);

data_rhu=importdata(filename_rhu);

data_ssd=importdata(filename_SSD);

data_win=importdata(filename_WIN);

data_prs=importdata(filename_prs);

totalsum=[];

for i=1:size(station,1)

sta=station(i);

syz_prec=find(data_prec(:,1)==sta);

syz_temp=find(data_temp(:,1)==sta);

syz_rhu=find(data_rhu(:,1)==sta);

syz_ssd=find(data_ssd(:,1)==sta);

syz_win=find(data_win(:,1)==sta);

syz_prs=find(data_prs(:,1)==sta);

if length(syz_prec)==length(syz_temp) && length(syz_temp)==length(syz_rhu) && length(syz_rhu)==length(syz_ssd) && length(syz_ssd)==length(syz_win) && length(syz_prs)==length(syz_win)

dates1=data_prec(syz_prec,1:7);

precz=data_prec(syz_prec,8:10);

tempz=data_temp(syz_temp,8:10);

rhuz=data_rhu(syz_rhu,8:9);

ssdz=data_ssd(syz_ssd,8);

winz=data_win(syz_win,8:9);

prsz=data_win(syz_prs,8);

total1=[dates1,precz,tempz,rhuz,ssdz,winz,prsz]; %1-7 8-10 11-13 14-15 16 17-18 ,19

totalsum=[totalsum;total1];

end

end

datasum=[datasum;totalsum];

end

%首先对异常数据进行处理

%台站海拔高度 +100000 当台站海拔高度为估测值时，在估测数据基础上加100000

dem=datasum(:,4);

dem(dem>100000)=dem(dem>100000)-100000;

datasum(:,4)=dem;

%各要素项 32766 数据缺测或无观测任务

datasum(datasum==32766)=NaN; %待插补数据

% 风速 +1000 当风速超过仪器测量上限时，在上限数据基础上加1000

winz=datasum(:,17:18);

winz(winz>1000)=winz(winz>1000)-1000;

datasum(:,17:18)=winz;

% 降水量 32700 表示降水"微量"

% 32XXX XXX为纯雾露霜

% 31XXX XXX为雨和雪的总量

% 30XXX XXX为雪量(仅包括雨夹雪，雪暴)

prec=datasum(:,8:10);

prec(prec==32700)=0.1;%相当于0.01mm了，原有的数据单位是0.1mm

prec(prec>32000 & prec<32999)=prec(prec>32000 & prec<32999)-32000;

prec(prec>31000 & prec<31999)=prec(prec>31000 & prec<31999)-31000;

prec(prec>30000 & prec<30999)=prec(prec>30000 & prec<30999)-30000;

datasum(:,8:10)=prec;

% -10000 实际温度(零下)超仪器下限刻度，在下限数据基础上减10000

temp=datasum(:,11:13);

temp(temp

datasum(:,11:13)=temp;

% 对缺测数据进行填补，采用线性内插法

total=[];

for i=1:size(station,1)

sta=station(i);

sy=find(datasum(:,1)==sta);

if length(sy)>=365

value=datasum(sy,:); %得到该站点全部的数据

cd=size(value,1);

% 依次进行插值处理

for k=8:19 %要素从第15列开始

element=value(:,k);

syz=find(~isnan(element)); %找到非NaN值

if length(syz)>340 & length(syz)

if isnan(element(1))

element(1)=element(syz(1)) % 如果第一个为NaN值，则用邻近值进行填补

end

if isnan(element(cd))

element(cd)=element(syz(end)) % 如果最后个为NaN值，则用邻近值进行填补

end

valuez=interp1(syz,element(syz),[1:cd])';%得到有效点

else

valuez=element;

end

value(:,k)=valuez;

end

end

total=[total;value];

end

filename=['F:\Day1951-2012yuanshi\xj\',int2str(year),'年新疆气象要素数据集.xlsx'];

xlswrite(filename,total)

end

结果如下

QQ图片20191228141625.png

注意：经度，纬度和高程的顺序在使用时请再次确认下，作者表示很懒，不想在去确认了

更多需求，请查看个人介绍