形成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文件
image.png
image.png
选择里面的world，在里面选择wgs84
image.png
结果如下图所示
image.png
（3）采用投影转换工具使得站点文件投影和dem投影保持一致。
image.png
image.png
红色部分设置成和dem一致
（4）在属性表中进行经纬度转换
第三部结果出来后能够得到下面红色部分
image.png
单击选中一列后右键计算几何，分别得到经纬度
image.png
转换后的结果如下：
image.png
除掉不在dem范围内的点
由于dem范围的缘故，可能导致部分站点不在dem的范围内，此时需要对这些站点就行删除，具体如下
首先打开dem范围，并将Nodata部分变成黄色，显示出真实的dem范围
image.png
image.png
通过上图可知最底下的两个点不在范围上，需要进行剔除，进入编辑状态即可进行删除，删除后才能进行下一步工作。
高程提取
高程作为协变量可以直接通过dem进行提取。
image.png
image.png
结果如下：
image.png
后续将得到每个站点每天的经纬度和高程信息。
更多需求，请查看个人介绍
