ANUSPLIN插值专题系列-插值脚本编写及后期处理
往期回顾
在前期的数据准备中我们获得了ANUSPLIN插值格式的excel文件，需要进一步通过SPSS转化为固定ASCII码格式的dat文件，并记录下dat文件各表头的格式，准备完成后就可以进行脚本的编制以及后期的数据处理。
因此本文的主要内容分为三个部分：
（1）Excel格式转化为dat格式；
（2）插值脚本的编制及修改说明；
（3）插值数据后期处理。
1. 打开spss后，加载进去excel数据；
2. 设置好对应的表头信息，如下图所示：
需要注意的是，站点的数据类型为字符型，设置好后，另存为固定ASCII码格式的：
在保存的时候注意spss会有个输出的日志文件，下图中标红的部分需要在脚本中用到。
设置好后，可以开始写插值的脚本了，脚本主要用到splina.exe和lapgrd.exe.两个程序。
Splina的脚本文件如下：
```
ssd2007
0 #插值的单位，0表示无单位
2 #自变量的个数，这里指经纬度
1 #协变量的个数，这里指dem，注意也可以将高程设置为自变量，不用协变量
0
0
1342316.1 6206316.3 0 1 ##前两数是对应dem文件中经度的范围（在arcgis中的source看），但要比dem中的范围要大一点才行
1649788.01 5937788.2 0 1 ##前两数是对应dem文件中经度的范围，但也要比dem中的范围要大一点才行
-163.422 7371.14 1 1 ##前两数dem文件的上下限
1000
0
3
365 #需要插值的个数，这是是对每天进行插值，故是365个
0
1
1
ssd2007.dat #插值的dat文件名
1000 #站点的个数，可以设置的比实际站点数量大些
5 #标签的宽度，我们的站点都是五位数的，所以是5
(a12,f12.4,f12.4,f12.2,365f12.1) #数据的格式，与上面spss的输出日志保持一致
ssd2007.res
ssd2007.opt
ssd2007.sur
ssd2007.lis
ssd2007.cov
(注意这里要留一行)
```
里面的#标注部分根据自己的数据来进行修改即可，```这个符号这里只是用来囊括插值代码的，在写入文件的时候不用写进去。将上述信息写在一个txt文件中，然后另存为一个cmd格式的文件。假设这里另存的文件名称为splina_ssd2007.cmd。Lapgrd的脚本文件如下：
```
ssd2007.sur
0 #表示输出所有的插值面
1
1
1
1342316.25165 6206316.25165 8000 #前两个是dem文件中的经度范围，8000：插值的分辨率
2
1649788.1524 5937788.1524 8000 #前两个是dem中的纬度范围，8000：插值的分辨率，最好和dem的分辨率一致
0
2
chinadem.txt #ascii码的dem文件名称
2
-9999.0
ssd2007_001.grd
ssd2007_002.grd
ssd2007_003.grd
ssd2007_004.grd
ssd2007_005.grd
ssd2007_006.grd
ssd2007_007.grd
ssd2007_008.grd
ssd2007_009.grd
.......
ssd2007_365.grd
#注意这里有个空行才行
```
同样将上述信息写在一个txt文件中，然后另存为一个cmd格式的文件。假设这里另存的文件名称为lapgrd_ssd2007.cmd。新建一个运行的脚本，命名为run.cmd
```
splina
lapgrd
```
将代码复制到run.cmd中，双击run.cmd就能够运行得到结果，注意：上述文件必须都放在同一个文件夹下。插值结果是grd格式的，需要在arcgis中批量转换为tif格式。
(1)基于python的grd格式批量转为tif格式
grd格式是传统的arcinfo的格式，不适合进一步的计算，需要将grd转换为tif格式，本文采用arcgis中自带的python来进行批量转换。
具体代码如下所示：
```
import arcpy
arcpy.env.workspace="D:\\chazhi\\" #存放数据的文件夹
a=arcpy.ListRasters("*","grd") #得到文件夹下所有的grd名称
for i in a:
arcpy.RasterToOtherFormat_conversion(i,"D:\\chazhi\\tif\\","TIFF")
```
将上述代码复制到自带的python上后运行即可。
(2)基于matlab的气象要素异常值处理
插值完后会发现气象要素如降水会存在着负值，负值是由于引进了协变量DEM造成的，可以认为出现的负值的地方降水为0的，尤其是在进入日尺度上插值时，当有n多个这样文件需要处理时，通过arcgis中的栅格计算器来一个个实现尤为费时，本文提供一个基于matlab的处理，同时也能够加入投影信息。
```
[a,R]=geotiffread('F:\项目\dem.tif');%先导入投影信息
info=geotiffinfo('F:\项目\dem.tif');
e=dir('*.tif');
for i=1:size(e,1)
data=importdata(e(i).name);
data(data<0)=0;
geotiffwrite(e(i).name,data,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end
```
在运行matlab时，需要将工作目录调整到刚刚转换成tif文件的目录下。
以上就是“ANUSPLIN插值专题系列”的全部内容了。接下来，对本专题系列进行一个简要的总结：
（1）气象数据的收集与整理，整理成ANUSPLIN插值所需格式，中间需要用到Excel、SPSS和Matlab，其中Matlab针对长时间序列复杂的气象数据整理；
（2）对气象站点经纬度的转化，将各个站点的经纬度转换成m为单位，需要用到AIbes投影的DEM文件；
（3）DEM文件转化为ASCII码；
（4）脚本文件的构造与修改；
（5）后期插值数据的进一步处理。
欢迎小伙伴们批评指正插值过程中出现的错误，后期会收集各种错误，针对各个错误进行详细的解答。
完
资源仅供学术交流使用，严禁商用！如有侵权，请联系小编微信：ylc107919293
推文期数：2018127
下期预告：考研资讯|近三年研究生就业率Top10
责任编辑：杨俊凤 尹礼唱 卢雅焱
推文审核：张天舒 梁龙武 骆丹云
总审核：学术无界顾问团
往期回顾：