基于ANUSPLIN的批量气象插值-从数据处理到最终结果（2）

2018-06-30 本文已影响343人画长空_yin

形成ausplin所需要的格式【1】

ausplin所需要的格式在（1）中已经提及到了，本文这次主要讲述如何得到这种格式

获取全部的站点信息

当我们有许多年的站点数据时，可能会发现每年的站点数量不太一样，可能2010年有820个站点，而到了2011年只有810个站点，这时我们首先需要得到全部的站点信息，为后续的匹配打下基础。

读取原始的txt文件的前三列数据，分别是站点，经度和纬度，对第一列求唯一值，既可以得到全部的站点信息了。代码如下：

%authour yinlichang3064@163.com

%得到全部txt的信息
filename='H:\Day1951-2012yuanshi\';%%根目录，自行修改
e=dir(); %得到文件夹的名称
termsum=[];
for i=5:10
    file=strcat(filename,e(i).name,'\');
    if i==5
        id='-13011-';
    end
    if i==6
        id='-10004-';
    end
    if i==7
        id='-13003-';
    end
    if i==8
        id='-14032-';
    end
    if i==9
        id='-12001-';
    end
    if i==10
        id='-11002-';
    end
    for m=2000:2015
        for n=1:12
            if n<10
                term=importdata(strcat(file,'SURF_CLI_CHN_MUL_DAY-',file(28:30),id,int2str(m),'0',int2str(n),'.TXT'));
            else
                term=importdata(strcat(file,'SURF_CLI_CHN_MUL_DAY-',file(28:30),id,int2str(m),int2str(n),'.TXT'));
            end
            term=term(:,1:3);
            termsum=[termsum;term];
        end
    end
end

%获取其中的唯一值
tersum_sta=termsum(:,1);%全部站点
tersum_lon=termsum(:,3);%经度
tersum_lat=termsum(:,2);%纬度
tersum_unique=unique(tersum_sta);
allsum=[];
for i=1:length(tersum_unique)
    data=tersum_unique(i);
    sy=find(tersum_sta==data);
    tersum_lon1=tersum_lon(sy(1));
    tersum_lat1=tersum_lat(sy(1));
    %由于上面经纬度是以度分为单位的，所以需要进一步转化为以度为单位的
    if  tersum_lon1>=10000
        tersum_lon1=int2str(tersum_lon1);
        tersum_lon2=str2num(tersum_lon1(1:3))+str2num(tersum_lon1(4:5))./60;
    else
        tersum_lon1=int2str(tersum_lon1);
        tersum_lon2=str2num(tersum_lon1(1:2))+str2num(tersum_lon1(3:4))./60;
    end
    
    if tersum_lat1>=10000
        tersum_lat1=int2str(tersum_lat1);
        tersum_lat2=str2num(tersum_lat1(1:3))+str2num(tersum_lat1(4:5))./60;
    else
        tersum_lat1=int2str(tersum_lat1);
        tersum_lat2=str2num(tersum_lat1(1:2))+str2num(tersum_lat1(3:4))./60;
    end
    all=[data,tersum_lon2,tersum_lat2];
    allsum=[allsum;all];%allsum即包含全部的站点信息且不重复
end
xlswrite('不重复的站点信息.xlsx'，allsum)

出来的结果如下图所示，总共有839个不重复的点

image.png

投影转换

由于ausplin需要的是以m为单位的经纬度，故我们需要在arcgis中进行投影转换得到m为单位的经纬度。首先需要准备以阿尔波斯投影的dem文件，后期该文件也是作为协变量参与插值的。
（1）首先加载dem文件
（2）在excel中手动添加表头，station，x 和y，然后通过add xy的方式加载到arcgis中，并重新输出为shapile文件