如何在 Python 中使用克里金法插入二维空间数据? [英] How to interpolate 2D spatial data with kriging in Python?

问题描述

我有一个空间二维域，比如 [0,1]×[0,1].在这个域中，有 6 个点可以观察到一些感兴趣的标量(例如，温度、机械应力、流体密度等).如何预测未观察点的兴趣量?换句话说，我如何在 Python 中插入空间数据?

例如，考虑以下二维域中点的坐标(输入)和相关数量(输出)的相应观测值.

将 numpy 导入为 np坐标 = np.array([[0.0,0.0],[0.5,0.0],[1.0,0.0],[0.0,1.0],[0.5,1.],[1.0,1.0]])观察 = np.array([1.0,0.5,0.75,-1.0,0.0,1.0])

X 和 Y 坐标可以通过以下方式提取:

x =坐标[:,0]y = 坐标[:,1]

以下脚本创建一个散点图，其中黄色(相应的蓝色)代表高(相应的)输出值.

将 matplotlib.pyplot 导入为 pltfig = plt.figure()plt.scatter(x, y, c=observations, cmap='viridis')plt.colorbar()plt.show()

我想使用克里金法来预测二维输入域内规则网格上感兴趣的标量.我如何在 Python 中执行此操作?

解决方案

在

我们看到元模型的预测等于观察到的输入点处的观察.

这个元模型是坐标的平滑函数:它的平滑度随着协方差核平滑度的增加而增加，平方指数协方差核恰好是平滑的.

I have a spatial 2D domain, say [0,1]×[0,1]. In this domain, there are 6 points where some scalar quantity of interest has been observed (e.g., temperature, mechanical stress, fluid density, etc.). How can I predict the quantity of interest at unobserved points? In other words, how may I interpolate spatial data in Python?

For example, consider the following coordinates for points in the 2D domain (inputs) and corresponding observations of the quantity of interest (outputs).

import numpy as np
coordinates = np.array([[0.0,0.0],[0.5,0.0],[1.0,0.0],[0.0,1.0],[0.5,1.],[1.0,1.0]])
observations = np.array([1.0,0.5,0.75,-1.0,0.0,1.0])

The X and Y coordinates can be extracted with:

x = coordinates[:,0]
y = coordinates[:,1]

The following script creates a scatter plot where yellow (resp. blue) represents high (resp. low) output values.

import matplotlib.pyplot as plt
fig = plt.figure()
plt.scatter(x, y, c=observations, cmap='viridis')
plt.colorbar()
plt.show()

I would like to use Kriging to predict the scalar quantity of interest on a regular grid within the 2D input domain. How can I do this in Python?

解决方案

In OpenTURNS, the KrigingAlgorithm class can estimate the hyperparameters of a Gaussian process model based on the known output values at specific input points. The getMetamodel method of KrigingAlgorithm, then, returns a function which interpolates the data.

First, we need to convert the Numpy arrays coordinates and observations to OpenTURNS Sample objects:

import openturns as ot
input_train = ot.Sample(coordinates)
output_train = ot.Sample(observations,1)

The array coordinates has shape (6, 2), so it is turned into a Sample of size 6 and dimension 2. The array observations has shape (6,), which is ambiguous: Is it going to be a Sample of size 6 and dimension 1, or a Sample of size 1 and dimension 6? To clarify this, we specify the dimension (1) in the call to the Sample class constructor.

In the following, we define a Gaussian process model with constant trend function and squared exponential covariance kernel:

inputDimension = 2
basis = ot.ConstantBasisFactory(inputDimension).build()
covariance_kernel = ot.SquaredExponential([1.]*inputDimension, [1.0])
algo = ot.KrigingAlgorithm(input_train, output_train,
                           covariance_kernel, basis)

We then fit the value of the trend and the parameters of the covariance kernel (amplitude parameter and scale parameters) and obtain a metamodel:

# Fit
algo.run()
result = algo.getResult()
krigingMetamodel = result.getMetaModel()

The resulting krigingMetamodel is a Function which takes a 2D Point as input and returns a 1D Point. It predicts the quantity of interest. To illustrate this, let us build the 2D domain [0,1]×[0,1] and discretize it with a regular grid:

# Create the 2D domain
myInterval = ot.Interval([0., 0.], [1., 1.])
# Define the number of interval in each direction of the box
nx = 20
ny = 10
myIndices = [nx-1, ny-1]
myMesher = ot.IntervalMesher(myIndices)
myMeshBox = myMesher.build(myInterval)

Using our krigingMetamodel to predict the values taken by the quantity of interest on this mesh can be done with the following statements. We first get the vertices of the mesh as a Sample, and then evaluate the predictions with a single call to the metamodel (there is no need for a for loop here):

# Predict
vertices = myMeshBox.getVertices()
predictions = krigingMetamodel(vertices)

In order to see the result with Matplotlib, we first have to create the data required by the pcolor function:

# Format for plot
X = np.array(vertices[:,0]).reshape((ny,nx))
Y = np.array(vertices[:,1]).reshape((ny,nx))
predictions_array = np.array(predictions).reshape((ny,nx))

The following script produces the plot:

# Plot
import matplotlib.pyplot as plt
fig = plt.figure()
plt.pcolor(X, Y, predictions_array)
plt.colorbar()
plt.show()

We see that the predictions of the metamodel are equal to the observations at the observed input points.

This metamodel is a smooth function of the coordinates: its smoothness increases with covariance kernel smoothness and squared exponential covariance kernels happen to be smooth.

这篇关于如何在 Python 中使用克里金法插入二维空间数据?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！