生态遥感的学习笔记

整理气象数据

2019-06-28  本文已影响10人  荔枝猪

需求

现有一个txt文件,其中包含了2003-2005年所有站点的日气象数据,需要从中提取并整理需要的温度数据。

代码

%% 读取txt文件中的站点,温度等数据
%By Yang,2019/5/12
%所有站点,所有数据都在一个txt中
%先利用fopen打开文件,赋予文件代号fid
fid = fopen('I:\climate\raw_climate\data2002-2005.txt');
%利用textscan读取txt里的数据
% %q为读取双引号括起来的
% 'Delimiter',',' 以,为分隔符
% 'headerlines', 1 不读第一行
data = textscan(fid,'%q%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,8};           %温度数据V12001
tempp5 = str2double(temp5);

%组合矩阵
temm = [tempp1,tempp2,tempp3,tempp4,tempp5];

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

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

%% 找到每个站点数据中的缺测值,并填补为NAN
%因为数据源有问题,缺测的值没有填充,被直接删除了,所以要先找到缺测的天数(位置)
perfect = tem02_05_d{1,1}(:,1:4);          %创建一个有完整序列站点(本数据为第一个站点)的前四列,站点号,年,月,日
for k = 1:length(tem02_05_d)               %删掉缺测站点所剩的站点数
    data = tem02_05_d{1,k};
    a = tem02_05_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
    tem02_05d{1,k} = perfect;
end

%% 处理特征值(本数据32766为缺测的特征值),将其转为NAN值
for i = 1:length(tem02_05d)  
    tem02_05d{1,i}(tem02_05d{1,i} == 32766) = NaN; 
end

%% 将2002-2005年数据中NAN值较多的站点删去
%统计每个站点元胞中的NaN值数量
for k = 1:length(tem02_05d)
    temp = tem02_05d{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<50    %准备把数据中缺少5%以上的删去
      B(1,j)=1;    %将超过5%以上NaN的数据位置标记为1
   end
end
[p,q]=find(B==0);  %[p,q]为空元胞所在行列号,即缺少5%以上的数据
tem02_05d(q)=[];   %去掉空元胞(即删去),这里考虑的是原始数据是1行n列的情况,如果是m行1列就改成A(p)=[];
station(q)=[];     %对应站点删掉

%% 对缺测值进行线性插值,利用fillmissing函数
for i = 1:length(tem02_05d)  
    tem02_05d{1,i}(:,5) = fillmissing(tem02_05d{1,i}(:,5),'linear');
end

%% 单位转换,本数据温度需除以10
for k = 1:length(tem02_05d)
    temp = tem02_05d{1,k}; 
    temp(:,5) = temp(:,5)/10;
    tem02_05d{1,k} = temp;
end
上一篇下一篇

猜你喜欢

热点阅读