how to get better Kriging result graphs in openturns?

Question

I performed spherical Kriging, but I can't seem to get good output graphs. The coordinates(x, and y) range from around around 51 latitude and around 6.5 as longitude my observations range from -70 to +10 here is my code :

import openturns as ot
import pandas as pd
# your input / output data can be easily formatted as samples for openturns


df = pd.read_csv("kreuzkerpenutm.csv")


inputdata = ot.Sample(df[['x','y']].values)
outputdata = ot.Sample(df[['z']].values)


dimension = 2  # dimension of your input (x,y)
basis = ot.ConstantBasisFactory(dimension).build()
covarianceModel = ot.SphericalModel(dimension)
    
algo = ot.KrigingAlgorithm(inputdata, outputdata, covarianceModel, basis)
algo.run()
result = algo.getResult()
metamodel = result.getMetaModel()

lower = [-10.0] * 2 # lower bound of the 2D window
upper = [50.0] * 2 # upper bound of the 2D window
graph = metamodel.draw(lower, upper)
graph.setBoundingBox(ot.Interval(lower, upper))
graph.add(ot.Cloud(inputdata)) # overlay a scatter plot of the observation points
graph.setTitle("Kriging metamodel")

# A View object allows us to interact with the underlying matplotlib figure 
from openturns.viewer import View
view = View(graph, legend_kw={'bbox_to_anchor':(1,1), 'loc':"upper left"})
view.getFigure().tight_layout()

here is my output:

kriging metamodel graph

I don't know why my graph won't show me my inputs aswell as my kriging results.

thanks for ideas and help

Answer 1

Sorry for the late answer.

Which version of openturns are you using? Probably you have an embedded transformation of (input) data, which makes the data range between (-3, 3) approximately (standard scaling). The kriging result should contains the transformation in such a case. With more recent openturns implementations, this feature has been removed.

Hope this can help. Cheers

Answer 2

If the input data is not scaled in [-1,1]^d, the kriging metamodel may have issues to identify the scale parameters using maximum likelihood optimization. In order to help for this, we may:

• provide a better starting point for the scale parameters of the covariance model (this is trick "A" below),
• set the bounds of the optimization algorithm so that the interval where the parameters are searched for correspond to the data at hand (this is trick "B" below).

This is what the following script does, using simulated data instead of a csv data file. In the script, I create the data using a g function which is scaled so that it produces results in the [-10, 70] range, as in your problem. Please look carefuly at the setScale() method which sets the initial value of the covariance model: this is the starting point of the optimization algorithm. Then look at the setOptimizationBounds() method, which sets the bounds of the optimization algorithm.

import openturns as ot


dimension = 2  # dimension of your input (x,y)
distribution = ot.ComposedDistribution([ot.Uniform(-10.0, 50.0)] * dimension)
inputdata = distribution.getSample(100)

g = ot.SymbolicFunction(["x", "y"], ["30 + 3.0 * sin(x / 10.0) * (y / 10.0) ^ 2"])

outputdata = g(inputdata)

basis = ot.ConstantBasisFactory(dimension).build()
covarianceModel = ot.SphericalModel(dimension)
covarianceModel.setScale(inputdata.getMax())  # Trick A
algo = ot.KrigingAlgorithm(inputdata, outputdata, covarianceModel, basis)
# Trick B, v2
x_range = inputdata.getMax() - inputdata.getMin()
scale_max_factor = 2.0  # Must be > 1, tune this to match your problem
scale_min_factor = 0.1  # Must be < 1, tune this to match your problem
maximum_scale_bounds = scale_max_factor * x_range
minimum_scale_bounds = scale_min_factor * x_range
scaleOptimizationBounds = ot.Interval(minimum_scale_bounds, maximum_scale_bounds)
algo.setOptimizationBounds(scaleOptimizationBounds)
algo.run()
result = algo.getResult()
metamodel = result.getMetaModel()
metamodel.setInputDescription(["x", "y"])
metamodel.setOutputDescription(["z"])

lower = [-10.0] * 2  # lower bound of the 2D window
upper = [50.0] * 2  # upper bound of the 2D window
graph = metamodel.draw(lower, upper)
graph.setBoundingBox(ot.Interval(lower, upper))
graph.add(ot.Cloud(inputdata))  # overlay a scatter plot of the observation points
graph.setTitle("Kriging metamodel")

# A View object allows us to interact with the underlying matplotlib figure
from openturns.viewer import View

view = View(graph, legend_kw={"bbox_to_anchor": (1, 1), "loc": "upper left"})
view.getFigure().tight_layout()

The previous script produces the following figure.

There are other ways to implement trick B. Here is one provided by J.Pelamatti:

# Trick B, v3
for d in range(X_train.getDimension()):
    dist = scipy.spatial.distance.pdist(X_train[:,d])
    scale_max_factor = 2.0  # Must be > 1, tune this to match your problem
    scale_min_factor = 0.1  # Must be < 1, tune this to match your problem
    maximum_scale_bounds = scale_max_factor * np.max(dist)
    minimum_scale_bounds = scale_min_factor * np.min(dist)

This topic is discussed in this particular thread in OT's forum.