如何在 Python 中使用克里金法插入测站数据?
Posted
技术标签:
【中文标题】如何在 Python 中使用克里金法插入测站数据?【英文标题】:How can I interpolate station data with Kriging in Python? 【发布时间】:2017-12-23 19:38:53 【问题描述】:浏览网页我发现在 Python 中使用克里金法的一些工具是 pyKriging 和 Gaussian Process Regression。但是,我无法让它们中的任何一个工作。第一个对我不起作用(甚至无法导入):
import pyKriging
File "~/python3.6/site-packages/pyKriging/krige.py", line 142
except Exception, err:
^
SyntaxError: invalid syntax
第二个我不明白如何使用它。我找不到一个简单的工作示例(例如 rroowwllaanndd answer 很好,但遗憾的是数据不再可供下载)
所以我的问题是,我如何使用克里金法插入我的数据?我有几个站数据保存在 numpy 数组中,如下所示:
2000 1 1 5.0
2000 1 2 3.4
2000 1 3 0.2
列是年 - 月 - 日 - 降水。我有几个这样的数据数组(st1、st2、st3)和另一个数组,其中包含每个站的 ID 和每个站所在的坐标(stid,所以站 1 位于经度 15.6865、纬度 62.6420 和等等)。
import numpy as np
st1 = np.array([[2000,1,1,5.0],[2000,1,2,3.4],[2000,1,3,0.2]])
st2 = np.array([[2000,1,1,8.2],[2000,1,2,2.5],[2000,1,3,0.0]])
st3 = np.array([[2000,1,1,np.nan],[2000,1,2,4.5],[2000,1,3,1.2]])
stid = np.array([[1,15.6865,62.6420],[2,15.7325,62.1254],[3,16.1035,61.1449]])
我需要的是每天的数组(或 3D 数组),其中包含每天在这样的网格中使用克里金法插值的所有站点的数据:
y = np.arange(61,63,0.125)
x = np.arange(14,17,0.125)
X,Y = np.meshgrid(x,y)
感谢任何帮助。
【问题讨论】:
【参考方案1】:很高兴找到有趣的文档、包等,克里金法通常被称为“高斯过程回归”。
在 python 中,有很多示例的一个很好的实现是著名的机器学习包scikit-learn 之一。它基于著名的 DACE matlab 实现。
documentation for Gaussian Process Regression 包括5 tutorials,以及list of available kernels。
使用您提供的数据,您只需执行以下操作即可将简单模型与您选择的内核相匹配:
import sklearn
gp = sklearn.gaussian_process.GaussianProcessRegressor(kernel=your_chosen_kernel)
gp.fit(X, y)
【讨论】:
声明gp.fit(X, y)
在这种情况下不合适,因为 X 是一个二维经度数组,而 y 是一个一维纬度数组。【参考方案2】:
使用OpenTURNS,KrigingAlgorithm
估计条件高斯过程的超参数。您需要的元模型的目的是将(经度,纬度)2D 点作为输入,将给定日期的降水作为输出。
第一步是准备数据。在以下脚本中,我创建了包含经度/纬度对的 coordinates_train
变量和包含降水量的 precipitation_train
变量。我使用了 2000 年 1 月 2 日的降水量,因为 2000 年 1 月 1 日在 3 站的数据缺失。
import openturns as ot
# Input points
coordinates_train = ot.Sample([[15.68,62.64],[15.73,62.12],[16.10,61.14]])
# Output points
precipitation_train = ot.Sample([[3.4],[2.5],[4.5]]) # At 2000/1/2
然后我们可以训练克里金法。为此,我使用常数基础(模型的趋势)和指数协方差模型。这应该是合适的,因为降水相对于站点的位置必须非常规律。
# Fit
inputDimension = 2
basis = ot.ConstantBasisFactory(inputDimension).build()
covarianceModel = ot.SquaredExponential([1.]*inputDimension, [1.0])
algo = ot.KrigingAlgorithm(coordinates_train, precipitation_train, covarianceModel, basis)
algo.run()
result = algo.getResult()
krigingMetamodel = result.getMetaModel()
然后我们可以使用元模型来预测未记录位置的降水量。由于krigingMetamodel
是一个函数,所以我只使用“()”运算符。
# Predict
coordinates = [15.70,62.53] # A new latitude/longitude pair
precipitation = krigingMetamodel(coordinates)
那么precipitation
是一个包含给定位置降水的一维点。这是预测的降水量。
>>> print(precipitation)
[3.46667]
您不妨通过将(经度、纬度、时间)作为输入来获得更一般的克里金法。在这种情况下,您所要做的就是向输入训练样本添加一个新维度,其中包含关联的时间,格式化为真实值。
【讨论】:
以上是关于如何在 Python 中使用克里金法插入测站数据?的主要内容,如果未能解决你的问题,请参考以下文章
使用 gstat 和 automap 包的克里金法 - 复制教程时的问题