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. ThegetMetamodel
method ofKrigingAlgorithm
, then, returns a function which interpolates the data.First, we need to convert the Numpy arrays
coordinates
andobservations
to OpenTURNSSample
objects:The array
coordinates
has shape(6, 2)
, so it is turned into aSample
of size 6 and dimension 2. The arrayobservations
has shape(6,)
, which is ambiguous: Is it going to be aSample
of size 6 and dimension 1, or aSample
of size 1 and dimension 6? To clarify this, we specify the dimension (1) in the call to theSample
class constructor.In the following, we define a Gaussian process model with constant trend function and squared exponential covariance kernel:
We then fit the value of the trend and the parameters of the covariance kernel (amplitude parameter and scale parameters) and obtain a metamodel:
The resulting
krigingMetamodel
is aFunction
which takes a 2DPoint
as input and returns a 1DPoint
. 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: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 thevertices
of the mesh as aSample
, and then evaluate the predictions with a single call to the metamodel (there is no need for afor
loop here):In order to see the result with Matplotlib, we first have to create the data required by the
pcolor
function:The following script produces the plot:
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.