基于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）