基于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 模块描述

1.3 ANUSPLIN 光滑样条函数模型选择

1.4 ANUSPLIN 插值模型选择

*注：对于日、年等数据，也可尝试选择以上模型。

二、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]编写说明

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]编写说明

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)