⽓象基于ANUSPLIN 的⽓象数据插值过程
  本⽂共分为3个部分，分别从 模型原理，⽓象数据处理，插值脚本编写 三个⽅向进⾏展开介绍。⽬标：  1、降低学习门槛，学会利⽤SPSS、Excel和ArcGIS⼯具实现ANUSPLIN插值的全流程。
  2、熟练使⽤SPSS（25.0）、Excel（2016）、ArcGIS（10.8）及ANUSPLIN（4.3）等软件。
  3、锻炼耐⼼、阅读⼒和理解⼒。本例以⽂字和表的形式介绍全部流程（重点部分加粗表⽰），学习时需要认真阅读和思考。
  当然，本⽂所介绍的⽓象数据整理⽅法并⾮为最简洁的⽅法，但是确实最为基础的⽅法。这对那些编程能⼒较弱（主要为⼊门程序员）和Excel能⼒⼀般的⼈员⽽⾔⾮常有⽤。例如：在进⾏⽓象站点数据提取时，Excel的VLOOKUP函数其实更加有效和简单，在将⽓象数据整理为所需样式时，Excel的数据透视表也能实现相应的功能。当然，这两步可能对基础较差的⼈来说较难理解。⽽本⽂所列举的步骤虽然不如上两步迅速快捷，但处理步骤更加连贯，环环相扣，容易理解，更易上⼿。当然，对于资深程序员来说，这些步骤统统可以⽤代码直接搞定，例如Python中的pandas库，如果程序写得好，即可在数秒内完成所有数据的清洗和转换；也可以写Python脚本来⽣成ANUSPLIN脚本；如果你的代码能⼒和数学能⼒都很强，也可以将ANUSPLIN算法直接⽤Python或其他编程语⾔直接复现，岂不是更易使⽤。
  本教程写于2018年本⼈研⼀期间，当时正处于各种学习积累的开始期，也不如现在⼀样对软件、编程有⼀定的了解。因此本教程即是我初学（不管是科研还是软件、编程能⼒）阶段探索过程的成果积累，也是让我慢慢到状态，⾛出迷茫的开始。  最后，希望所有读者有所收获，也希望能与各位同僚⼤神交流经验，共同学习，共同进步。
⽬录
0 ANUSPLIN 模型的使⽤
0.1 模型算法
  ANUSPLIN 采⽤局部薄盘光滑样条法。其理论统计模型表达⽅法如下：
  是位于空间  点的因变量；为  维样条独⽴变量， 是要估算的关于 的未知光滑函数； 为  维独⽴协变量；  为 的  维系数；  为具有期望值为 0 且⽅差为 的⾃变量随机误差； 是作为权重的已知局部相对变异系数， 为误差协⽅差，在所有数据点上为常数，但通常未知。
  若 ，则模型简化为简化为薄盘光滑样条原型。函数 和系数 采⽤最⼩⼆乘法估计来确定：
  是函数的粗糙度测度函数(样条次数)，定义为函数 的  阶偏导，是正的光滑参数，通常由⼴义交叉验证 GCV 的最⼩化来确定。
  其他更具体的参数计算参见 ANUSPLIN 说明⽂件。
0.2 ANUSPLIN 模块描述
表 0.2 ANUSPLIN 模块描述
程序模块
描述
z =i f (x )+i b y +T i e  i =i
1,2,…,N
z i  i  x i  d  f  x i  y i  p  b  y i  p  e i  w σi 2 wi  σ2 p =0 f  b ((z −i =1∑
N
i f (x )−i b y )/w )+T i i 2ρJ (f )
m  J (f )m  f (x )i  f  m  ρ
SPLINA适⽤于站点数<2000 的任意个独⽴变量或多个协变量的薄盘样条函数。数据平滑度由 GCV 或 GML 决定。SPLINB功能与 SPLINA 类似，站点数>2000 数据利⽤ SELNOT 进⾏节点选择，最多 10000 个站点，2000 个节点。SELNOT为 SPLINB 添加初始节点
ADDNOT为 SPLINB 添加数据节点
DELNOT为 SPLINB 删除数据节点
GCVGML 计算每个表⾯的 GCV 或 GML，以及平滑参数范围内所有表⾯的平均GCV 或 GML。可⽤于样条优化参数。将 GCV 或 GML 值写⼊⽂件，以
供检查和绘图。
LAPPNT计算预测值或贝叶斯标准误差估计的点⽂件
LAPGRD⽣成拟合曲⾯或贝叶斯标准误差曲⾯
程序模块描述
*注：本例对SPLINA、LAPGRD模块的详细参数进⾏解释。
0.3 ANUSPLIN 光滑样条函数模型选择
表 0.3 18 个薄盘光滑样条函数模型详细列表
模型序号变量（协变量）样条次数模 型 缩写含义1经度，纬度2BVTPS2双变量薄盘光滑样条函数
2经度，纬度3BVTPS3双变量薄盘光滑样条函数
3经度，纬度4BVTPS4双变量薄盘光滑样条函数
4经度，纬度(⾼程)2TVPTPS2三变量局部薄盘光滑样条函数
5经度，纬度(⾼程)3TVPTPS3三变量局部薄盘光滑样条函数
6经度，纬度(⾼程)4TVPTPS4三变量局部薄盘光滑样条函数
7经度，纬度,⾼程(m)2TVTPS2三变量薄盘光滑样条函数
8经度，纬度,⾼程(m)3TVTPS3三变量薄盘光滑样条函数
9经度，纬度,⾼程(m)4TVTPS4三变量薄盘光滑样条函数
10经度，纬度,⾼程(km)2TVTPS2三变量薄盘光滑样条函数
11经度，纬度,⾼程(km)3TVTPS3三变量薄盘光滑样条函数
12经度，纬度,⾼程(km)4TVTPS4三变量薄盘光滑样条函数
13经度，纬度,⾼程(dm)2TVTPS2三变量薄盘光滑样条函数
0.4 ANUSPLIN 插值模型选择
表 0.4 不同⽓象要素空间插值的最佳模型参考
⽓象要素模型独⽴变量独⽴协变量数据转换⽅式样条次数⽉平均最⾼⽓温TVPTPS经度，纬度⾼程2,3⽉平均最低⽓温TVPTPS经度，纬度⾼程2,3⽉平均风速TVPTPS经度，纬度⾼程2⽉降⾬量BVTPS经度，纬度-平⽅根转换2⽔汽压TVPTPS经度，纬度⾼程2⽉⽇照时数TVPTPS经度，纬度⽉温度范围2,3
净辐射SRAD+TVPTS
经度，纬度年辐射，地表覆盖，地形等
2,3蒸发⽫蒸发
QVPTPS
经度，纬度
净辐射,⽔⽓压差，风速
2,3,4
⽓象要素模型独⽴变量独⽴协变量数据转换⽅式
样条次数*注：对于⽇、年等数据，也可尝试选择以上模型。
1 数据整理
1.1 提取所需站点信息
  ⾸先需要获取所拥有的⽓象数据中的⽓象站点信息以及研究区域的*.shp ⽂件，参考坐标系尽量相同。利⽤ ARCGIS（版本为 10.5）或其他⽅法进⾏站点筛选，筛选原则为研究区内所有站点和研究区域边界外垂直于边界⽅向 1-2 个站点的站点数的总和。参考⽅法如下：  1）导⼊数据。在 ARCGIS 中，点击<；添加数据>，导⼊包含站点经纬度数据的 .xls ⽂件和对应的sheet。*注：ARCGIS不能很好地打开.xlsx ⽂件，建议将包含站点经纬度数据的.xlsx ⽂件另存为.xls ⽂件。
  2）显⽰数据。右键点击 <；图层> 下导⼊的⽂件数据，选择 <；显⽰ XY 数据> ， <X 字段> 选择经度列， <Y 字段> 选择纬度列，<Z 字段> 可不进⾏设置，<；输⼊坐标的坐标系> 选择 WGS 1984 。
  3）将数据⽂件导出为 .shp ⽮量图层⽂件。⽅法：[右键数据-数据-导出数据]。
  4） 导⼊研究区域的.shp ⽮量图层⽂件，并对其进⾏编辑。利⽤ <；编辑器> ⼯具，按照筛选原则，制作⼀个包含所有所需站点的.shp ⽮量图层⽂件。
  5）裁剪所需要的点数据。以制作的包含所有所需站点的.shp ⽮量图层⽂件为掩膜， 裁剪导出的⽓象站点.shp ⽂件，得到所需的⽓象站点。⽅法为：[地图处理-裁剪]，输⼊要素为站点导出的.shp ⽮量图层⽂件，裁剪要素为制作的包含所有所需站点的.shp ⽮量图层⽂件。
  6）导出数据。打开 [地图处理-ArcToolBox-转换⼯具-表转 Excel] ，将裁剪后得到的⽓象站点数据⽂件导出为 Excel ⽂件，并删除或替换其中的冗余数据。
1.2 原始⽓象数据整理
1.2.1 ⽓象数据整理格式
表 1.2.1a ANUSPLIN 插值数据格式⽰例
区站号经度纬度⾼程2018.1.1
2018.1.2
2018.1.3
2018.1.4
50136122.5252.97438.5253550246124.7252.35361.9147650349124.451.67501.5046650353126.6351.73173.9066750425120.1850.25581.4536650434
121.68
50.48
732.6
5
3
7
7
*注：此处经纬度坐标系为 WGS-1984 地理坐标系，单位为度（°），也可为 albers 投影坐标系，单位为⽶ （m），或其他坐标系。经纬度⾼程为变量信息，可进⾏取舍，如果还有其他变量或协变量，可在经纬度⾼程后依次添加相应变量数据。
表 1.2.1b 原始⽓象数据格式⽰例
区站号
纬度
经度
⾼程
年
⽉
⽇
平均⽓温最⾼⽓温最低⽓温平均⽓温质量控
迷你世界二郎神图片制码⽇最⾼⽓温质量控
制码⽇最低⽓温质量控
制码5013652581223143301999121-205-165-2420005013652581223143301999122-293-175-36100050136
5258
12231
4330
1999
12
3
-270
-173
-324
∗∗∗∗∗∗∗∗∗∗∗
5013652581223143301999124-248-136-3140005013652581223143301999125-230-117-29600050136
5258
12231
4330
1999
12
6
-230
oracle连接驱动名
language is the bridge 演讲-115
-306
区站号纬度经度⾼程年⽉⽇平均⽓温最⾼⽓温最低⽓温平均⽓温质量控
制码⽇最⾼⽓温质量控
制码⽇最低⽓温质量控
制码1.2.2 原始数据合并
  原始⽓象站点数据为*.TXT 格式⽂本，若⽂本较多，需要将其合并从⽽简化步骤。以国家⽓象数据共享⽹的⽓象站点数据为例，通过<SURF_CLI_CHN_MUL_DAY_FORMAT.doc>《数据集内容和⽂件组织⽅式》⽂件信息，出所需要的⽓象数据对应的.txt ⽂件。例如⽇照时数数据，以<SURF_CLI_CHN_MUL_DAY-SSD-14032-YYYYMM.TXT>的⽅式命名。  合并步骤如下：
  1）将所需⽓象数据拷贝⾄空⽂件夹，在此⽂件夹下新建⼀个.txt ⽂件。  2）打开新建的.txt ⽂件，编辑以下内容：
type  *.txt > <⽂件名>.txt
  其中<⽂件名>为合并后⽂件的⽂件名。
  3）保存⽂件，将⽂件扩展名 txt 改为 bat，双击运⾏。⽂件夹下所有.txt ⽂件会合并到<⽂件名>.txt ⽂件。*注：原始.txt ⽂件末尾必须有空⾏，否则会出现合并乱⾏。
表 1.2.2 “⽇照时数”数据组织⽅式⽰例
序号中⽂名数据类型单位
1区站号Number(5)2纬度Number(5)（度、分）3经度Number(6)（度、分）4观测场海拔⾼度
oracle对比两个表字段Number(7)0.1⽶5年Number(5)年6⽉Number(3)⽉7⽇Number(3)⽇8⽇照时数Number(7)0.1⼩时9
⽇照时数质量控制码
Number(2)
*数据解析⽅式为：每⾏的字符串长度为（5+5+6+7+5+3+3+7+2），与数据类型列（）内标记的数据长度对应且相等。其中1-5位为区站号，6-10位为纬度……依此类推。
1.2.3 数据重构
  数据重构利⽤ SPSS 软件（版本为 25.0）。步骤如下：
  1）打开 SPSS，利⽤ [⽂件-导⼊数据-⽂本数据] 导⼊合并的⽓象数据⽂件。
*注：若⽓象数据含有特定特征值 32766 或其他特殊值⽽导致相邻数据之间⽆空格符区分，则可能导致 SPSS ⽆法正常读取这些数据。这⾥，可以利⽤ Excel 来解决这些问题，以⽇照时数数据为例：
  1. 依据《数据集内容和⽂件组织⽅式》⽂件，获取对应的序号、中⽂名、数据类型、单位等信息， 提取数据类型 Number 后的数据宽度值信息。
  2. 打开 Excel，导⼊合并后的⽓象数据⽂件，导⼊数据选项中选择固定宽度，依据 Number 后的数据宽度值划分数据，完成后另存数据，利⽤ SPSS 导⼊另存后数据，也可在此处直接进⾏< 2）>的步骤后再进⾏ SPSS 导⼊数据。
∗∗∗∗∗
  2）检查数据信息，删除⽆效数据。此处可保留经纬度⾼程数据，若经纬度⾼程数据相同站点出现不同值，建议剔除经纬度⾼程，保留⽇期和所需的⽓象数据即可。气象python零基础入门教程
*建议去除经纬度⾼程信息，利于<1.2.4 站点筛选>步骤。
  3）数据重构。⽅法为：[数据-重构] (可不进⾏保存) 。
1.2.4 站点筛选
  若重构后的数据完全为所研究的数据，则忽略此步骤。若只利⽤其中部分站点的数据，则利⽤ SPSS 进⾏以下步骤：  1）获得所需的站点数据，包括经纬度⾼程信息。  2）导⼊重构后的⽓象数据和站点数据⽂件。
  3）⽂件合并。将⽓象数据合并到⽓象站点数据，步骤为 [数据-合并⽂件-添加变量] ， 打开⽓象数据的数据集，选择基于键值的⼀对⼀合并，键值选择区站号，点击确定。将合并后的⽂件另存为 Excel ⽂件。
  4）删除多余站点数据。打开导出的 Excel ⽂件，筛选出没有经纬度⾼程信息的站点， 删除站点及⽓象数据，并保存。
1.2.5 缺失值插补
  在 SPSS 中进⾏缺失值插补。利⽤同⼀站点不同⽇期数据，对缺失数据通过线性内插的⽅法进⾏数据补全。⽅法如下：  1）利⽤ SPSS 或 Excel 对筛选过的站点数据进⾏转置。保证每⼀列为同⼀站点不同⽇期的数据。  2）将特定特征值 32766 或其他特殊值替换为空值，并删除经纬度⾼程信息。span元素
  3）打开 [转换-替换缺失值] ⼯具，导⼊第⼀组要插值的数据，⽅法选择 <；线性内插> ， 之后点击 <；变化量>，最后，选择其余所有需要缺失值插补的数据组导⼊，点击确定后完成差值。
  4）删除插值前的数据，保留插值后的数据，重新转置数据并导出为 Excel ⽂件。
1.2.6 更换单位并改变⼩数位数
  依据《数据集内容和⽂件组织⽅式》，各个⽓象数据的单位并不为 1。例如降⽔的单位为 0.1mm，即，如果降⽔数据为 135，则真实降⽔量应为 135 * 0.1 = 13.5mm，另外，缺失值插补后⼩数位数会增加。若单位和⼩数位数不影响后续计算和分析，则忽略此步，否则，则可利⽤进⾏以下步骤进⾏单位换算和改变⼩数位数：
  1）打开缺失值插补后 Excel ⽂件，新建⼀个sheet，复制第⼀⾏⽇期和第⼀列区站号⾄新sheet 相应位置。  2）从新sheet 中对应的原sheet 第⼀个数据的位置开始，输⼊以下公式：
=ROUND (Sheet1!B2*0.1,1)
  其中 Sheet1!B2 为原 sheet 第⼀个数据的位置，0.1 为数据单位，1 为需要保留的⼩数位数。*注意：在 Excel 中单元格设置中设置⼩数位数，只能改变显⽰的⼩数位数，源数据⼩数位数并未改变。
  3）将公式填充到新 sheet 中所有对应数据的位置，完成单位换算和⼩数位数的设置。建议将新sheet 中的数据，以⽂本的格式复制到另⼀个⽆数据的sheet 或新 Excel ⽂件中，便于下⼀步使⽤。
1.2.7 重新添加站点经纬度和⾼程信息
  如果需要转换坐标系，可以利⽤<*1.3 站点坐标转换>的⽅法进⾏经纬度坐标转换。  由于前步骤去除了站点经纬度⾼程信息，需要在数据整理最后步骤中重新进⾏添加。
  添加过程如<1.2.4 站点筛选>过程中的⽂件合并，基本思路为，将处理完成的⽓象数据合并到所需站点数据之后，由于站点相同，所以合并完成的⽂件便为包含站点经纬度⾼程信息和⽓象数据的⽂件。也可将经纬度⾼程信息直接粘贴到数据⽂件的相应位置，同样可以得到<；表 1.2>所⽰例的 Anusplin 插值数据格式。最后导出数据⽂件，建议为 Excel ⽂件。
1.2.8 导出⽓象数据⽂件
  ANUSPLIN 要求输⼊数据⽂件为以 ASCII 码保存的包含站点经纬度和⾼程信息的站点数据（.dat），可以使⽤ SPSS 进⾏数据⽂件导出。
  1）SPSS 导⼊包含完整数据的 Excel ⽂件。将数据第⼀⾏读取为变量名称。*注：导出的.dat 不含变量名称。
∗∗

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。