简书大神Anusplin气象插值的逐步实现（GIS小白重现此过程的一些步骤总结，请大家多指正）

上学期师兄介绍了一个简书大神画长空_yin分享的ANUSPLIN气象插值法，昨天捣鼓了一晚上，试了好几次，也重新看了大神的博客和网友的提问与答疑，终于在今天上午试成功了，现在想记录下自己操作的过程，供大家参考。
先放上简书领路大神的博客供大家学习，包括了前期数据处理等好几个详细步骤，我本人主要参考和学习的是（2）和（5）：
基于ANUSPLIN的批量气象插值-从数据处理到最终结果（1） - 简书 https://www.jianshu.com/p/07ccc981404d
基于ANUSPLIN的批量气象插值-从数据处理到最终结果（2） - 简书 https://www.jianshu.com/p/c6979cf3f619
基于ANUSPLIN的批量气象插值-从数据处理到最终结果（3） - 简书 https://www.jianshu.com/p/ccb07139ce56
基于ANUSPLIN的批量气象插值-从数据处理到最终结果（4) - 简书 https://www.jianshu.com/p/2869eaa5c208
基于ANUSPLIN的批量气象插值-从数据处理到最终结果（5） - 简书 https://www.jianshu.com/p/f3ed50f0ffea
基于ANUSPLIN的批量气象插值-从数据处理到最终结果（6） - 简书 https://www.jianshu.com/p/df78282f25d6

首先，说明下本研究的需求：需要将北京、天津和廊坊3市的41个气象站点 1999-2018年逐月的降雨数据进行空间插值。
接下来就是插值，一共分为5步。

第一步：准备所需数据
包括气象数据和高程数据

(1)气象站点信息及气象数据（本研究主要涉及的降水数据）
这是我初步整理的41个气象站点240个月的逐月数据，其中，num_station是台站编号，x和y分别经纬度（单位为度）

先将表格（要是97-2003版本的表格）加载进GIS，以坐标xy显示，并导出shp格式的文件、加载进GIS。然后定义坐标（本研究是京津廊，所以选的UTM北半球50，可以根据自己的研究需求设定，但要与接下来设置的dem的投影保持一致）
加载表格数据

转化为shp文件加载进来，定义坐标


打开shp文件的属性表，将x和y值重新计算成以米为单位的经纬度（点击x和y所在列，右键Calculate Geometry）

这是已重新计算的经度x，可以看到经度x变成了以“米”为单位的数据，但纬度y还是以“度”为单位，接下来就把y也变为以米为单位的数即可。

(2)高程数据
加载已切好的1km网格的高程dem数据（京津廊地区的dem），并将投影转化为与上述气象矢量一致的北半球50

然后，提取dem中的数据到气象shp文件，并生成新的气象shp文件（即比以前多一列属性数据，就是dem）

自此，所需的两类数据就准备好了，得到气象站点的编号num_station、经纬度x和y（单位是米）、高程dem（在dem的栅格文件里提取的数），并将除了站点编号外的其他数据的数据格式设置为“保留一位小数”（方便后面SPSS导出ASCII码数据时数据格式统一），如下表所示。

第二步：转化数据
还是包括气象数据和高程数据

（1）气象数据“20200425.dat”文件生成
将上面的excel表加载到spss软件里（操作：打开数据-类型选excel）

在此窗口对数据进行下检查和设置，主要是将台站编码num_station的类型设为“字符串”，并检查下时间序列的气象数据的数据格式是否统一，比如我刚才设置都是保留一位小数，导入spss后这里显示的小数位数就是1，此外，最好把度量标准除了台站编码外都设置为“度量”类型。
然后，将数据另存为ASCII格式备用（即20200425.dat）。此外，在关闭SPSS之前，需将操作日志导出成txt格式，后面编脚本时需要用到台站编码、经纬度、高程和气象数据的数据格式等信息，如A11、F11.1等

（2）高程数据“dem_jjl_new.txt”文件生成
用GIS下面这个工具，将高程数据转化为ASCII码的文件（即dem_jjl_new.txt）

第三步：编写脚本
编写脚本，一共要编写两个splina需要运行的脚本ssdpres和lapgrd需要运行的脚本ssdprel（名字可以根据自己的需求进行变更）

（1）splina需要运行的脚本设置如下（图片转自大神）：

需要注意最后面不是留一行，需要留8行，这是一位网友试出来的，我师兄和我的插值操作也是留了8行才成功的。

下面是我自己的脚本，注意#和#后面的文字需要删除，那个只是大神帮助大家说明下怎么修改。可以看到有些要素是根据我自己的需求来修改的，如文件名称ssdpre、经纬度的范围、气象数据的dat文件、数据的格式等。

（2）lapgrd需要运行的脚本设置如下（图片转自大神）：

我自己的脚本如下，注意#和#后面的文字需要删除，那个只是大神帮助大家说明下怎么修改。可以看到有些要素是根据我自己的需求来修改的，如文件名称ssdpre、经纬度的范围、高程数据的txt文件等，以及我的需求是240个月的数据，所以要把所有需要的grd都列出来（从ssdpre1999_01一直到ssdpre2018_12）。

将上述的ssdpres和ssdprel的txt文件的后缀改为.cmd文件，得到ssdpres.cmd和ssdprel.cmd

（3）总运行脚本设置如下（图片转自大神）：
新建一个名为run.txt，写入如下代码，并将后缀改为cmd，即run.cmd：
我的代码如下：
splina<ssdpres.cmd>ssdpres.log
lapgrd<ssdprel.cmd>ssdprel.log

第四步：运行脚本

（1）把如下文件放进一个文件夹：
气象数据 20200425.dat
高程数据 dem_jjl_new.txt
总运行脚本 run.cmd
splina运行脚本 ssdpres.cmd
lapgrd运行脚本 ssdprel.cmd
splina运行程序 splina.exe(网上可以下载）
lapgrd运行程序lapgrd.exe(网上可以下载）

（2）双击run.cmd
如果运行成功后，可以看到你的文件夹里多出来了这些文件，其中grd后缀的就是插值的月降水数据。

如果没有正常运行，可看下简书大神博客下面的答疑，以及生成的两个日志log，里面可以看到代码运行到哪一步

第五步：批量将GRD格式转化为TIFF格式
用“栅格转其他”这个工具进行
Input rasters 加入需要转化的批量的grd文件，output workspace选择导出的tiff保存的文件夹或数据库，此外，别忘了在environment setting中设置投影（与本研究保持一致北半球50），就ok了