基于ANUSPLIN的批量气象插值-从数据处理到最终结果(3)

在(2)中我们已经得到了想要的全部不重复的站点信息,在(3)部分我将展示如何将(1)中的数据结合(2)来转变成我们ausplin所需要的插值格式。代码如下所示:

a=xlsread('D:\日降水插值\中国气象站点shp\china-station1.xlsx'); %(2)中得到的站点信息
filename='H:\Day1951-2012yuanshi\插值格式\中国';
station=a(:,1);x=a(:,2);y=a(:,3);dem=a(:,4);%站点,经度,纬度,高程;
e=dir(fullfile('*.xlsx'));% (1)中存放数据的目录
i=1;while i<=length(e)
    a=xlsread(e(i).name);
    st=a(:,1);year=a(:,2);month=a(:,3);day=a(:,4);termtemmax=a(:,5);
    sy=[];
    for j=1:length(station)
        sy1=find(st==station(j));
        sy=[sy;sy1];
    end
    st=st(sy);year=year(sy);month=month(sy);day=day(sy);termtemmax=termtemmax(sy);
    sy=find(termtemmax>=32000 & termtemmax<32766);
    termtemmax(sy)=0.1; %降水中的异常值》=32000 & <32766的是极小值,用0.1替代
    sta=unique(st);
    sta=sort(sta);%升序排列
    year1=year(1);
    if mod(year1,4)==0
        kk=366;
    else
        kk=365;
    end
    term1=zeros(length(sta),kk);
    for m=1:length(sta)
        sy1=find(st==sta(m));
        if length(sy1)==kk 
        term1(m,:)=termtemmax(sy1)';
        else
        term1(m,:)=zeros(1,kk)+32766;
        end
    end
    term1(find(term1==32766))=NaN;
    zz1=[];
    for n=1:length(sta)
        st1=sta(n);
        sy=find(station==st1);
        zz=[x(sy),y(sy),dem(sy)];
        zz1=[zz1;zz];
    end
    term2=[sta,zz1,term1];
    name1=e(i).name;
    name1=strcat(filename,name1);
    xlswrite(name1,term2);
    i=i+1;
end

上述代码将气象数据中的异常值32766以空值来表示,结果输出到excel时,异常值的部分是缺失的。

结果如下所示:


image.png

当一年中的异常值部分缺失较多时,则可以直接去掉该站点,而当一年的异常值缺失部分较少时,则可以通过线性插值法或邻近站点插值法来进行填补,具体见下期(4),上面代码在matlab上进行

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

推荐阅读更多精彩内容