使用Matlab处理站点气候数据

%% 读取txt文件中的站点，温度，降水,湿度等数据
% By Yang，2019/5/12
% 所有站点，所有数据都在一个txt中
% 先利用fopen打开文件，赋予文件代号fid
fid = fopen('I:\climate\raw_climate\气温降水湿度1961-2018\2018new.csv');
% 利用textscan读取txt里的数据
% %q为读取双引号括起来的数据,共九列
% 'Delimiter',',' 以,为分隔符
% 'headerlines', 1 不读第一行
data = textscan(fid,'%q%q%q%q%q%q%q%q','Delimiter',',','headerlines', 1);
fclose(fid);

temp1 = data{1,1}; %站点号V01000
tempp1 = str2double(temp1); %转换字符串为数值
temp2 = data{1,2}; %年V04001
tempp2 = str2double(temp2);
temp3 = data{1,3}; %月V04002
tempp3 = str2double(temp3);
temp4 = data{1,4}; %日V04003
tempp4 = str2double(temp4);
temp5 = data{1,6}; %降水数据V13201，温度数据V12001,平均相对湿度数据V13003...
tempp5 = str2double(temp5);

begin_year=2018; %读取数据的起始年份和终止年份
end_year=2018;
%组合矩阵
temm = [tempp1,tempp2,tempp3,tempp4,tempp5];

%% 按站点分成不同的矩阵
station = unique(temm(:,1) ); %提取不重复站点号
%不同站点的数据分为不同元胞，排为一行
for i = 1:length(station)
index = find(temm(:,1)==station(i));%利用find找到对应站点数据的行号
data_d{1,i} = temm(index,:);
end

%% 记录数据缺少太多的站点位置，并删掉
%创建一个完整的日历持续时间数组；caldays表示间隔天数
t = datetime(begin_year,01,01):caldays(1):datetime(end_year,12,31);
[y,m,d] = ymd(t); %返回t中日期时间值的年、月和天数
d_max = length(d); %总天数(数据长度)

B = zeros(1,length(data_d)); %创建相同数据长度的0矩阵
for j = 1:length(data_d)
number = length(data_d{1,j}); %每个站点元胞内的数据长度
if number>round(d_max*0.9) %准备把数据中缺少5%以上的删去
B(1,j)=1; %将没有缺少的位置标记为1
end
end
[p,q]=find(B==0); %[p,q]为空元胞所在行列号，即缺少5%以上的数据
data_d(q)=[]; %去掉空元胞(即删去),这里考虑的是原始数据是1行n列的情况，如果是m行1列就改成A(p)=[];
station(q)=[]; %对应的站点号也删掉

%% 找到每个站点数据中的缺测值，并填补为NAN
%因为数据源有问题，缺测的值没有填充，被直接删除了，所以要先找到缺测的天数（位置）
perfect = [y',y',m',d']; %组合为一个矩阵,第一列随便填充，后三列年，月，日
for k = 1:length(data_d) %删掉缺测站点所剩的站点数
data = data_d{1,k};
a = data_d{1,k}(1,1); %a为对应站点的站点号
perfect(:,1) = a; %第一列替换为对应站点的站点号，年月日不变
perfect(:,5) = NaN; %第五列填充为NaN
for i =1:length(data)
year = data(i,2);
month = data(i,3);
day = data(i,4);
%利用find找到年月日都正确的站点数据的行号
ind = find(perfect(:,2)==year&perfect(:,3)==month&perfect(:,4)==day);
perfect(ind,5) = data(i,5); %缺测的数据会没有填充即为NaN
end
data_d2{1,k} = perfect;
end

%% 处理特征值
%降雨32700为微量，将其设为0；32744，32766等为缺测，将其转为NAN值
for i = 1:length(data_d2)
data_d2{1,i}(data_d2{1,i} == 999999) = NaN; %温度等数据则不需考虑该步
end

%% 将数据中NAN值较多的站点删去
%统计每个站点元胞中的NaN值数量
for k = 1:length(data_d2)
temp = data_d2{1,k};
index = isnan(temp(:,5)); %利用isnan来确定是否为NaN，是为1，否为0
ID = find(index == 1); %利用find找到NaN的行号
ID_temp{1,k} = temp(ID,:); %ID_temp为每个NAN值的数量
end

% 将数据中NaN值太多的站点删掉
[m,n] = size(ID_temp);
B = zeros(1,n);
for j=1:n
number = length(ID_temp{1,j});
if number<round(d_max*0.05) %准备把数据中缺少5%以上的删去
B(1,j)=1; %将超过5%以上NaN的数据位置标记为1
end
end
[p,q]=find(B==0); %[p,q]为空元胞所在行列号，即缺少5%以上的数据
data_d2(q)=[]; %去掉空元胞（即删去）,这里考虑的是原始数据是1行n列的情况，如果是m行1列就改成A(p)=[];
station(q)=[]; %对应站点删掉

%% 对缺测值进行线性插值，利用fillmissing函数
for i = 1:length(data_d2)
data_d2{1,i}(:,5) = fillmissing(data_d2{1,i}(:,5),'linear');
end
%如果是降雨等数据，注意去掉负值
for i = 1:length(data_d2)
data_d2{1,i}(data_d2{1,i} < 0) = 0;
end

%% 单位转换
for k = 1:length(data_d2)
temp = data_d2{1,k};
temp(:,5) = temp(:,5)/100; %本数据降水需除以10
data_d2{1,k} = temp;
end