基于ANUSPLIN的气象数据插值
这篇文章是对ANUSPLIN这个插值工具进行简单的介绍，项目demo可以参考：
https://github.com/leeyang1991/ANUSPLIN
这个项目已经把从数据转换到脚本运行等一系列工作都用python实现了。
至于ANUSPLIN中的一些细节和参数说明，参考自：
https://blog.csdn.net/weixin_42155937/article/details/120673718
零、ANUSPLIN 介绍
空间插值：将空间上离散点的测量数据转换为连续的数据曲面的过程
ArcGIS软件中，有趋势面法、自然邻域法、样条函数法、克里金法等插值方法，但不能引入协变量
协变量：温度插值引入高程数据、降水量引入海岸线
ANUSPLIN 基于薄盘样条函数理论 ,引入多个影响因子作为协变量进行气象要素空间插值 ,大大提高插值精度 ,且 能同时进行多个表面的空间插值, 对时间序列的气象要素更加适合。
——《专用气候数据空间插值软件ANUSPLIN 及其应用 》
将 ANUSPLIN 的局部薄盘样条插值结 果分别与反向距离权重法和普通克吕格法的插值结果进行对比 ， 显示 ANUSPLIN 软件的插值误差最小 。 结果同样表明 ， 适当增加站点数量、 提高DEM精度可进一 步提高 ANUSPLIN软件 的插值精度。
——《基于 ANUSPLIN 软件的逐日气象要素插值方法应用与评估 》
一、ANUSPLIN 模型说明
1.1 模型算法
ANUSPLIN 采用局部薄盘光滑样条法。其理论统计模型表达方法如下：
$z_i = f(x_i)+b^T y_i + e_i i = 1 , 2 , … , N $
\(z_i\)是位于空间 \(i\)点的因变量； \(x_i\)为 \(d\) 维样条独立变量，\(f\)是要估算的关于 \(x_i\)的未知光滑函数； \(yi\)为 \(p\)维独立协变量； \(b\)为 \(y_i\) 的\(p\)维系数； \(e_i\) 为具有期望值为 0 且方差为 $w_i σ_2 $的自变量随机误差； \(w_i\)是作为权重的已知局部相对变异系数，\(σ2\) 为误差协方差，在所有数据点上为常数，但通常未知。
若 \(p=0\)，则模型简化为简化为薄盘光滑样条原型。函数\(f\)和系数\(b\)采用最小二乘法估计来确定：
\(\displaystyle\sum_{i=1}^{N}((z_i-f(x_i)-b^Ty_i)/w_i)^2 +ρJ_m(f)\)
\(J_m(f)\) 是函数 \(f(x_i)\) 的粗糙度测度函数(样条次数),定义为函数 \(f\) 的m阶偏导， ρ是正的光滑参数，通常由广义交叉验证 GCV 的最小化来确定。
其他更具体的参数计算参见 ANUSPLIN 说明文件。
1.2 模块描述
| 程序模块 | 描述 |
|---|---|
| SPLINA | 适用于站点数<2000 的任意个独立变量或多个协变量的薄盘样条函数。数据平滑度由 GCV 或 GML 决定。 |
| SPLINB | 功能与 SPLINA 类似，站点数>2000 数据利用 SELNOT 进行节点选择,最多 10000 个站点，2000 个节点。 |
| SELNOT | 为 SPLINB 添加初始节点 |
| ADDNOT | 为 SPLINB 添加数据节点 |
| DELNO | 为 SPLINB 删除数据节点 |
| DELNOT | 为 SPLINB 删除数据节点 |
| GCVGML | 计算每个表面的 GCV 或 GML，以及平滑参数范围内所有表面的平均GCV 或 GML。可用于样条优化参数。将 GCV 或 GML 值写入文件，以供检查和绘图。 |
| LAPPNT | 计算预测值或贝叶斯标准误差估计的点文件 |
| LAPGRD | 生成拟合曲面或贝叶斯标准误差曲面 |
1.3 ANUSPLIN 光滑样条函数模型选择
| 模型序号 | 变量（协变量） | 样条次数 | 模 型 缩写 | 含义 |
|---|---|---|---|---|
| 1 | 经度,纬度 | 2 | BVTPS2 | 双变量薄盘光滑样条函数 |
| 2 | 经度,纬度 | 3 | BVTPS3 | 双变量薄盘光滑样条函数 |
| 3 | 经度,纬度 | 4 | BVTPS4 | 双变量薄盘光滑样条函数 |
| 4 | 经度,纬度(高程) | 2 | TVPTPS2 | 三变量局部薄盘光滑样条函数 |
| 5 | 经度,纬度(高程) | 3 | TVPTPS3 | 三变量局部薄盘光滑样条函数 |
| 6 | 经度,纬度(高程) | 4 | TVPTPS4 | 三变量局部薄盘光滑样条函数 |
| 7 | 经度,纬度,高程(m) | 2 | TVTPS2 | 三变量薄盘光滑样条函数 |
| 8 | 经度,纬度,高程(m) | 3 | TVTPS3 | 三变量薄盘光滑样条函数 |
| 9 | 经度,纬度,高程(m) | 4 | TVTPS4 | 三变量薄盘光滑样条函数 |
| 10 | 经度,纬度,高程(km) | 2 | TVTPS2 | 三变量薄盘光滑样条函数 |
| 11 | 经度,纬度,高程(km) | 3 | TVTPS3 | 三变量薄盘光滑样条函数 |
| 12 | 经度,纬度,高程(km) | 4 | TVTPS4 | 三变量薄盘光滑样条函数 |
| 13 | 经度,纬度,高程(dm) | 2 | TVTPS2 | 三变量薄盘光滑样条函数 |
1.4 ANUSPLIN 插值模型选择
| 气象要素 | 模型 | 独立变量 | 独立协变量 | 数据转换方式 | 样条次数 |
|---|---|---|---|---|---|
| 月平均最高气温 | TVPTPS | 经度,纬度 | 高程 | 2,3 | |
| 月平均最低气温 | TVPTPS | 经度,纬度 | 高程 | 2,3 | |
| 月平均风速 | TVPTPS | 经度,纬度 | 高程 | 2 | |
| 月降雨量 | BVTPS | 经度,纬度 | - | 平方根转换 | 2 |
| 水汽压 | TVPTPS | 经度,纬度 | 高程 | 2 | |
| 月日照时数 | TVPTPS | 经度,纬度 | 月温度范围 | 2,3 | |
| 净辐射 | SRAD+TVPTS | 经度,纬度 | 年辐射,地表覆盖,地形等 | 2,3 | |
| 蒸发皿蒸发 | QVPTPS | 经度,纬度 | 净辐射,水气压差,风速 | 2,3,4 |
降水的空间插值由于数据范围大以及不确定性强，所以采用平方根变换来降低数据的值域范围，最后再运用普通双变量的薄盘样条函数。
*注：对于日、年等数据，也可尝试选择以上模型。
二、ANUSPLIN 使用方法
2.1 splina.cmd脚本
rainsp
1
2
1
0
0
140 150 0 5
-45 -39 0 5
2 5000 1 1
1000.0
2
2
12
0
1
1
rainn.dat
300
6
(a6,f8.3,f8.3,f8.1/12f9.2)
04-P-04-splina.res
04-P-04-splina.opt
04-P-04-splina.sur
04-P-04-splina.lis
04-P-04-splina.cov
注：每个参数为一行。
[splina.cmd]编写说明
| 参数（例） | 说明 | 解释 | |
| rainsp | 表面文件名 | ||
| 1 | 插值数据的单位 | 0-未定义 | |
| 1-米 | |||
| 2-英尺 | |||
| 3-千米 | |||
| 4-英里 | |||
| 5-度 | |||
| 6-弧度 | |||
| 7-毫米 | |||
| 8-兆焦耳 | |||
| 2 | 变量 | 独立样条变量 | 0-10 |
| 1 | 独立协变量 | 0-7 | |
| 0 | 表面样条变量 | 0-7 | |
| 0 | 表面协变量 | 0-7 | |
| 140 150 0 5 | 第1个变量（经度）的参数 | 含义： 前两个数字：经纬度和高程的边界信息。 第三个数字：表示经纬度和高程（变量）是否进行转换和转换方法。 第四个数字：表示单位。 如果有，第五个数字代表边距。 | ①经纬度和高程的边界（大于DEM中的范围） [140 150]-经度的西东边界 [-45 -39]-纬度南北边界 [2 5000]-高程的下上边界 ②经纬度 (变量)的变换方法 0-不进行变换 1-X/A 2-X*A 3-A*LOG(X + B) 4-(X/B)**A 5-A*EXP(X/B) 6-A*TANH(X/B) 7-取异向角值 8-取异向角系数 ③单位定义可参考[插值数据的单位] ④当第3个变量为高程时，可将单位m变换为km，可提高拟合效果。 |
| -45 -39 0 5 | 第2个变量（纬度）的参数 | ||
| 2 5000 1 1 | 第3个变量（协变量-高程）的参数 | ||
| 1000.0 | 变换系数A | A=1000.0 | |
| 2 | 因变量转换 | 0-不转换 | |
| 1-转换为自然对数 | |||
| 2-转换为平方根 | |||
| 2 | 样条次数 | ≥2 | |
| 12 | 输出表面数 | ≥1 | |
| 0 | 误差权重选择 | 0-所有数据点采用统一的权重 | |
| 1-所有的表面采用统一的权重 | |||
| 12(输出表面数)-为每个表面采用不同的权重 | |||
| 1 | 优化参数指标(通常为1) | 0-每个表面采用相同的平滑参数 | |
| 1-每个表面采用相同的平滑算法 | |||
| 2-每个表面采用不同的平滑算法 | |||
| 1 | 平滑方法选择 (通常为1) | 0-采用表面平滑参数-参数大小 | |
| 1-最小化GCV法 | |||
| 2-使用所提供的误差标准偏差估计最小化真实均方误差(MSE) | |||
| 3-采用固定自由度-固定值大小 | |||
| 4-最小化GML法 | |||
| rainn.dat | 插值数据文件 | ||
| 300 | 数据点个数 | ≥数据点真实值 | |
| 6 | 站点标签字符数 | ||
| (a6,f8.3,f8.3,f8.1/12f9.2) | 数据格式 | a6-标签类型为字符串，6字符。 f8.3-变量数据，3位小数，8个字符。 /-换行。 12f9.2-12个9个字符保留两位小数的数据。 若有： x-空格。 48x, f8.1,40x-数据为8个字符1位小数，除去前6(48/8)列数据和后5(40/8)列数据，只保留第7列数据。 （注意：定义的数据数量与站点标签、变量个数与需要插值的数据数量总和尽量一致。例如根据 <表1.2.1 ANUSPLIN插值数据格式示例> 的示例数据，我们可以定义 （a5, 2f8.2, f8.1, 4f8.0），表示1个5字符的站点标签，2个8字符2位小数的变量（经度、纬度），1个8字符1位小数的协变量（高程），4个8字符无小数的数据。） | |
| 若有，此处可以为导入导出数据节点文件和坏数据标志文件。若无，则忽略，不需空行。（脚本示例中为 忽略） | |||
| 04-P-04-splina.res | 输出残差文件 | 不需要可不填，需空行。 | |
| 04-P-04-splina.opt | 光滑参数文件 | 不需要可不填，需空行。 | |
| 04-P-04-splina.sur | 表面文件 | LAPGRD的输入数据 | |
| 04-P-04-splina.lis | 列表文件 | 观测数据、预测数据、误差数据。不需要可不填，需空行。 | |
| 04-P-04-splina.cov | 拟合表面系数的误差协方差文件 | 不需要可不填，需空行。 | |
| (建议空两行作为结尾) | |||
2.2 lapgrd.cmd脚本
04-P-04-splina.sur
1
1
1
04-P-04-splina.cov
2
1
1
144.5 148.5 0.25
2
-44.0 -40.0 0.25
0
2
tas4.txt
2
-99.0
rain1.grd
(f8.2)
2
-99.0
cov_rain1.grd
(f8.2)
[lapgrd.cmd]编写说明
| 参数（例） | 说明 | 解释 | ||
| 04-P-04-splina.sur | 输入表面名 | |||
| 1 | 表面个数 | 0-全部输出 | ||
| 1-12([splina.cmd]中的输出表面数)-输出指定个数的表面，与最后输出的*.grd个数一致。 | ||||
| 1 | 表面类型计算 | 0-只统计摘要信息 | ||
| 1-拟合表面值 | ||||
| 1 | 表面值转换，不一定需要。依[splina.cmd]中的因变量转换确定。没有则忽略，不需空行。 | 0-不进行转换 | ||
| 1-应用转换 | ||||
| 04-P-04-splina.cov | 误差协方差文件名，没有则忽略，不需空行。 | |||
| 2 | 误差协方差类型。没有则忽略，不需空行。 | 0-只计算平均表面值的标准误差 | ||
| 1-模型标准误差 | ||||
| 2-预测标准误差 | ||||
| 3-95%模型置信区间 | ||||
| 4-95%预测置信区间 | ||||
| 最大标准误差，可不填，需空行。 | ||||
| 1 | 栅格位置 | 0-栅格边角 | ||
| 1-栅格中心 | ||||
| 1 | 第1个变量的索引(经度) | 数据中第1数字类型的数据列 | ||
| 144.5 148.5 0.25 | 边界与分辨率，小于[splina.cmd]中纬度的边界范围，与DEM的范围一致。 | [144.5 148.5]-西东边界 0.25-分辨率 | ||
| 2 | 第2个变量的索引(纬度) | 数据中第2数字类型的数据列 | ||
| -44.0 -40.0 0.25 | 边界与分辨率。小于[splina.cmd]中纬度的边界范围，与DEM的范围一致。 | [-44.0 -40.0]-南北边界 0.25-分辨率 | ||
| 0 | 掩膜方式 | 0-无掩膜 | ||
| 1-通用掩膜 | ||||
| 2-ARCGIS掩膜 | ||||
| 3-Idrisi掩膜 | ||||
| 若有掩膜，此处可以为掩膜文件名。若没有，则忽略，不需空行。（脚本示例中为忽略） | ||||
| 2 | 独立协变量数据格式 | 0-固定常数 | ||
| 1-通用栅格格式 | ||||
| 2-ArcGIS格式 | ||||
| 3-Idrisi格式 | ||||
| tas4.txt | 独立协变量数据文件名。若为常量，此处为常量值。 | 此处为ArcGIS导出以ASCII码保存的纯文本dem数据。 | ||
| 2 | 输出文件格式 | 0-X,Y,Z网格文件 | ||
| 1-以行的形式保存通用栅格文件 | ||||
| 2-ARCGIS栅格文件 | ||||
| 3- Idrisi影像文件 | ||||
| -99.0 | 输出文件无效值定义 | 设置值与DEM中无效值相同。 | ||
| rain1.grd | 输出文件名 | 输出文件个数与[表面个数]一致。每个文件名占一行。 | ||
| (1f8.2) | 输出文件数量（1，为1时可忽略）和数据格式（f8.2） | 空白表示以二进制形式输出。 *注意：如果输出文件无效值设置过大，建议增加输出文件字符宽度。 | ||
| 2 | 输出误差表面格式 | 同[输出文件格式]。 | ||
| -99.0 | 输出误差文件无效值定义 | 同[输出文件无效值定义]。 | ||
| cov_rain1.grd | 输出误差文件名 | 同[输出文件名]。 | ||
| (1f8.2) | 输出误差文件数量（1，为1时可忽略）和数据格式（f8.2） | 同[输出文件数量和数据格式]。 | ||
| (建议空两行作为结尾) | ||||
2.3 grd转tif文件
import arcpy
from arcpy.sa import *
import sys, string, os
idx = 'Avst'
dir = r'E:\projects\weather\ANUSPLIN\mapdata\\'
path = dir + idx
files = os.listdir(path)
for f in files:
if os.path.splitext(f)[1] == '.grd':
Input_raster_file = path + os.sep + f
Raster_Format = "TIFF"
Output_Workspace = path
basename = os.path.splitext(f)[0]
Output_raster = Output_Workspace + os.sep + basename + '.tif'
arcpy.RasterToOtherFormat_conversion(Input_raster_file, Output_Workspace, Raster_Format)
print(Output_raster)
os.remove(Input_raster_file)
2.4 负值问题
有些气象要素如降水、辐射等插值完后会可能会存在着负值，不符合逻辑，可以将出现的负值的地方设为0，方法可参考下面文章。
https://www.jianshu.com/p/e7ae1e367f8a
import arcpy
from arcpy.sa import *
arcpy.env.workspace = "I:\\climate\\vpd\\tif_02_05" #工作空间（文件夹）
output_path = "I:\\climate\\vpd\\tif05\\" #输出地址
rasterlist = arcpy.ListRasters("*","tif") #遍历文件夹下所有的tif格式的文件
for raster in rasterlist:
out = output_path + raster
outCon = Con(Raster(raster) < 0, 0,raster)
outCon.save(out)
