How to get a better krigging result

[ChironeX1976]

I have been studying the tutorials on the gstools website (and others).
I posted this question in the pykrige section too, but i think it does not belong there.

I succeeded in producing a plot of the speed in a cross section of a river based on krigging.
It is still wrong in certain points.
How can i force the method to get a value in a certain point?
Should i take another krigging method?
Please, help. It's for the fish in there :-)
ugly heatmap, before krigging:

nicer heatmap, after krigging:

my new code:

import plotly.graph_objs as go
from plotly.subplots import make_subplots
import numpy as np
import gstools as gs
from matplotlib import pyplot as plt
# data preparation: x y z values are a cross section of the speed of riverflow
riverbed = np.array([[0.15, 0.], [0.3, -0.2], [0.5, -0.25], [0.8, -0.32], [1.1, -0.42], [1.3, -0.42], [1.5, -0.42],
                     [1.7, -0.33], [1.8, -0.28], [1.9, -0.07], [2., 0.]])
xyz = np.array([[ 2.,0.,0.],[ 1.9, 0.07,0.002], [1.9,0.21,-0.002],[ 1.9,0.28,-0.006],[ 1.9,0.35,0.],[ 1.8,0.066,0.052],
       [ 1.8,0.198,0.02 ],[ 1.8,0.264,0.012],[ 1.8,0.33,0.],[ 1.7,0.07,0.138],[ 1.7,0.21,0.11 ],[ 1.7,0.28,0.097],
       [ 1.7,0.35,0.], [ 1.6,0.082,0.121], [ 1.6,0.246,0.097],[ 1.6,0.328,0.071],[ 1.6,0.41,0.],[ 1.5,0.084,0.046],
       [ 1.5,0.252,0.035],[ 1.5,0.336,0.073],[ 1.5,0.42,0.],[ 1.3, 0.084,0.005],[ 1.3,0.252,0.018],[ 1.3,0.336,0.021],
       [ 1.3,0.42,0.],[ 1.1,0.084,0.018],[ 1.1,0.252,-0.002],[ 1.1,0.336,-0.004],[ 1.1,0.42,0.],[ 0.8,0.064,0.012],
       [ 0.8,0.192,0.007], [ 0.8,0.256,0.002],[ 0.8,0.32,0. ],[ 0.5,0.05,0.004], [0.5,0.15,-0.002], [ 0.5,0.2,-0.001],
       [ 0.5,0.25,0.], [ 0.3,0.04,0.002],[ 0.3,0.16,-0.004], [ 0.3,0.2,0.], [ 0.15,0.,0.]])
x=xyz[:,0]
y=-xyz[:,1]
z=xyz[:,2]
#plot ugly heatmap
# ------------------
fig = make_subplots()
fig.add_trace(go.Heatmap(x=x,y=y, z=z, type='heatmap', colorscale='bluered', zsmooth=False))
fig.add_trace(go.Scatter(x=x,y=y, mode='markers+text',text=z, name='speed', showlegend=True,
                         marker=dict(size=20, color = z, colorscale = 'bluered', showscale = False )))
fig.add_trace(go.Scatter(x=riverbed[:,0],y=riverbed[:,1],mode='lines'))
fig.show()
#########################
# 2. Estimation of the variogram
#########################
# Choose i_afstand and i_binstap based
# on visual interpretation of data (plot)
# fast testplot: to estimate nominal distance from graph
plt.scatter(x, y)
plt.xlabel('m from side')
plt.ylabel('depth')
plt.show()
# based on estimation of the parameters i_afstand and i_binstap (by iteration)
# you should get a slope in the datapoints getting to a plateau:
# left: correlated distances, right: uncorrelated distances
# start by estimating, plot and change till the variogram is right
# estimation of distance of correlated values
i_afstand = 0.55
# make pairs of numbers going from zero to a distance with discrete steps
i_binsstap = 0.1
bins = np.arange(start=0, stop=i_afstand, step=i_binsstap)
bin_center, gamma = gs.vario_estimate(pos=(x, y),
                                      field=z,
                                      bin_edges=bins)
# fast testplot: the variogram
plt.scatter(bin_center, gamma)
plt.xlabel('distance between pairs')
plt.ylabel('gamma')
plt.show()
# variogram is estimated? next step.
# ##############################################################################################
# # 3. Fit variogram with model
#   model based on "gs.Stable" function
#   Wackernagel, H. “Multivariate geostatistics”, Springer, Berlin, Heidelberg (2003)
# ##############################################################################################
#  create model object with 2 dimensions, and standard parameters, based on Stable-function
model = gs.Stable(dim=2)
# fit the variogram with a stable model.
fit_parameters, pcov, r2_score = model.fit_variogram(bin_center, gamma,
                                                     nugget=False,
                                                     return_r2=True)
# fast testplot: fitted curve with data
plt.figure()
ax = model.plot(x_max=max(bin_center))
ax.scatter(bin_center, gamma)
plt.xlabel('afstand tussen punten')
plt.ylabel('gamma')
plt.show()
# variogram is fitted with model? Next step.
##################################################
# # 4. Kriging based on model
###################################################
# grid definition, based on size of river
gridx = np.arange(0.15, 2, 0.2) # start, einde, stap
gridy = np.arange(-0.0, -0.45, -0.2)
OK = gs.krige.Ordinary(model, [x, y], z, exact=True)
OK.structured([gridx, gridy])
# fast testplot for kriged data "river is upside down"
plt.figure()
ax = OK.plot()
ax2 = plt.scatter(x,y)
plt.show()
# # plot nicer heatmap
# # ------------------
fig = make_subplots()
fig.add_trace(go.Heatmap(x=gridx,y=gridy, z=OK.all_fields[0].transpose(), type='heatmap', colorscale='bluered', zsmooth='fast'))
# rivierbed - line
fig.add_trace(go.Scatter(x=x,y=y, mode='markers+text',text=z, name='speed', showlegend=True,
                         marker=dict(size=20, color = z, colorscale = 'bluered', showscale = False )))
fig.add_trace(go.Scatter(x=riverbed[:,0],y=riverbed[:,1],mode='lines'))
fig.show()

[LSchueler]

Hi ChironeX1976,

for the fish, I'll jump right into it, before I go to bed ;-)

Kriging methods should (normally) always match the given observation points. Without having had a detailed look into your code, my first guess would be that the colour scales of the heatmap and of the scatter plots do not match, because the krige field and the obs. points have different minima and maxima. I have never used plotly before, but in matplotlib it would be something like setting vmin and vmax.

Based on the colour scale, the red area looks really dangerous for the fish ;-)

But in all seriousness, your data set looks really interesting and if you are interested, we could create a small kriging tutorial with your data.

[ChironeX1976]

Hi LSchueler,

thanks for your response. Great idea, that krigging tutorial in 2D with real data.

You are right about the scales. I updated my code, just before where i make the "nicer heatmap".

But i have a bigger problem: i think i withdrew the wrong kriged values from the OK-object.
Which values should i use?

# ... see code in previous post...
# fast testplot for kriged data "river is upside down"
plt.figure()
ax = OK.plot()
ax2 = plt.scatter(x,y)
plt.show()
# get z-value from OK
krigz = OK.all_fields[0].transpose()
# # plot nicer heatmap
# # ------------------------------------
fig = make_subplots()
# krigvalues (heatmap)
fig.add_trace(go.Heatmap(x=gridx,y=gridy, z=krigz, type='heatmap', zsmooth='fast',coloraxis = "coloraxis"))
# original values
fig.add_trace(go.Scatter(x=x,y=y, mode='markers+text',text=z, name='speed', showlegend=True, marker=dict(size=20, color = z, coloraxis="coloraxis")))
# rivierbed - line
fig.add_trace(go.Scatter(x=riverbed[:,0],y=riverbed[:,1],mode='lines'))
# make uniform colorscale for both datasets based on min and max values
cmin = min(z.min(), krigz.min())
cmax = max(z.max(), krigz.max())
fig.update_layout(coloraxis = dict(colorscale='bluered', cmin=cmin, cmax = cmax))
fig.show()

[LSchueler]

Hi again,
sorry, I was quite busy yesterday.
I think we really did not make it very clear how to access the kriged field as a data structure. After the kriging, it is stored in OK.field, which maybe looks a bit nicer than OK.all_fields[0]

I just made a quick plot of the resulting kriged field from your code (with a different colormap) and I guess your discrepancies stem from the low resolution of the field. This gets hidden by the colour interpolation of plotly.

You have created a structured grid on which you did the kriging. Instead, you can use an unstructured grid and hand over a list of positions, exactly like the positions x and y of your data. This way, you can precisely control at which positions you want the kriged field to be calculated. This way you could check that the kriged field matches your observations.
But you could also simply increase the resolution of the structured grid, see the 2nd figure.

I did this by add following to your code:

x_no = 100
y_no = 100
gridx = np.linspace(0.15, 2.0, x_no)
gridy = np.linspace(-0.0, -0.45, y_no)
OK.structured((interp_gridx, interp_gridy))
krigz = OK.field.T
# SNIP
fig.add_trace(go.Heatmap(x=gridx, y=gridy, z=krigz, type='heatmap', colorscale='bluered', zsmooth='fast'))

You can then simply increase your resolution further by increasing x_no and y_no.
And if you are still skeptical (looking at the 0.073 and 0.071 values near the riverbed), you should probably decrease the size of the scatter plot dots, as they have a uniform colour and hide the colour gradient of the kriged field beneath them.

I hope that helps!

I think we should update some of the tutorials to show how to actually work with the kriged data.

Cheers,
Lennart

[ChironeX1976]

Thanks again! You solved the problem. I share the complete, working code below.
I worked out the two options you came up with: 1. structured grid, 2. random grid.
The random grid makes more sense to me because the measurements were also done on a random grid.

structured grid

from plotly.subplots import make_subplots
import numpy as np
import gstools as gs
from matplotlib import pyplot as plt
# ##############################################################################
# 1. data preparation: x y z values are a cross section of the speed of riverflow
# ##############################################################################
riverbed = np.array([[0.15, 0.], [0.3, -0.2], [0.5, -0.25], [0.8, -0.32], [1.1, -0.42], [1.3, -0.42], [1.5, -0.42],
                     [1.6, -0.41],[1.7, -0.35], [1.8, -0.33], [1.9, -0.35], [2., 0.]])

xyz = np.array([[2.,0.,0.],[1.9,0.07,0.002],[1.9,0.21,-0.002],[1.9,0.28,-0.006],[1.9,0.35,0.],[1.8,0.066,0.052],
                [1.8,0.198,0.02],[1.8,0.264,0.012],[1.8,0.33,0.],[1.7,0.07,0.138],[1.7,0.21,0.11],[1.7,0.28,0.097],
                [1.7,0.35,0.],[1.6,0.082,0.121],[1.6,0.246,0.097],[1.6,0.328,0.071],[1.6,0.41,0.],[1.5,0.084,0.046],
                [1.5,0.252,0.035],[1.5,0.336,0.073],[1.5,0.42,0.],[1.3,0.084,0.005],[1.3,0.252,0.018],[1.3,0.336,0.021],
                [1.3,0.42,0.],[1.1,0.084,0.018],[1.1,0.252,-0.002],[1.1,0.336,-0.004],[1.1,0.42,0.],[0.8,0.064,0.012],
                [0.8,0.192,0.007],[0.8,0.256,0.002],[0.8,0.32,0.],[0.5,0.05,0.004],[0.5,0.15,-0.002],[0.5,0.2,-0.001],
                [0.5,0.25,0.],[0.3,0.04,0.002],[0.3,0.16,-0.004],[0.3,0.2,0.],[0.15,0.,0.]])
x=xyz[:,0]
y=-xyz[:,1]
z=xyz[:,2]


# #########################################
# # 2. Estimation of the variogram
# Choose i_afstand and i_binstap based
# on visual interpretation of data (plot)
# #########################################
#fast testplot: to estimate nominal distance from graph
plt.scatter(x, y)
plt.xlabel('m from side')
plt.ylabel('depth')
plt.show()

# based on estimation of the parameters i_afstand and i_binstap (by iteration)
# you should get a slope in the datapoints getting to a plateau:
# left: correlated distances, right: uncorrelated distances
# start by estimating, plot and change till the variogram is right
# estimation of distance of correlated values
i_afstand = 0.55
# make pairs of numbers going from zero to a distance with discrete steps
i_binsstap = 0.1
bins = np.arange(start=0, stop=i_afstand, step=i_binsstap)
bin_center, gamma = gs.vario_estimate(pos=(x, y),
                                      field=z,
                                      bin_edges=bins)
# fast testplot: the variogram
plt.scatter(bin_center, gamma)
plt.xlabel('distance between pairs')
plt.ylabel('gamma')
plt.show()

# variogram is estimated? next step.
# ##############################################################################################
# # 3. Fit variogram with model
#   model based on "gs.Stable" function
#   Wackernagel, H. “Multivariate geostatistics”, Springer, Berlin, Heidelberg (2003)
# ##############################################################################################

#  create model object with 2 dimensions, and standard parameters, based on Stable-function
model = gs.Stable(dim=2)
# fit the variogram with a stable model.
fit_parameters, pcov, r2_score = model.fit_variogram(bin_center, gamma,
                                                     nugget=False,
                                                     return_r2=True)
# fast testplot: fitted curve with data
plt.figure()
ax = model.plot(x_max=max(bin_center))
ax.scatter(bin_center, gamma)
plt.xlabel('afstand tussen punten')
plt.ylabel('gamma')
plt.show()

# variogram is fitted with model? Next step.
##################################################
# # 4. Kriging based on model
###################################################
# very small, structured, grid definition, based on size of river
x_no = 100 # interpolation steps
y_no = 100
gridx = np.linspace(0.15, 2.0, x_no)
gridy = np.linspace(-0.0, -0.45, y_no)

OK = gs.krige.Ordinary(model, [x, y], z, exact=True)
OK.structured([gridx, gridy])

# fast testplot for kriged data "river may be upside down"
plt.figure()
ax = OK.plot()
ax2 = plt.scatter(x,y)
plt.show()

# get z-value from OK
krigz = OK.field.transpose()

# # plot heatmap with plotly
# # --------------------------
fig = make_subplots()
# krigvalues (heatmap)
fig.add_trace(go.Heatmap(x=gridx,y=gridy, z=krigz,
                         type='heatmap', zsmooth='fast',coloraxis = "coloraxis"))
# original values
fig.add_trace(go.Scatter(x=x,y=y, mode='markers+text',text=z, name='speed', showlegend=True,
                         marker=dict(size=20, color = z, coloraxis="coloraxis")))
# riverbed - line
fig.add_trace(go.Scatter(x=riverbed[:,0],y=riverbed[:,1],mode='lines'))

# make uniform colorscale for both datasets based on min and max values
cmin = min(z.min(), krigz.min())
cmax = max(z.max(), krigz.max())
fig.update_layout(coloraxis = dict(colorscale='bluered', cmin=cmin, cmax = cmax))
fig.show()

random grid

# see other code above
# #################################################
# # 4. Kriging based on model
###################################################
# make very small, UNstructured, grid definition, based on size of river
# ----------------------------------------------------------------------
# initalise randomizer
seed = gs.random.MasterRNG(19970221)
rng = np.random.RandomState(seed())
# The measurements were unstructured, in fractions of a meter (float).
# The gstools randomizer generates whole positive numbers (integers).
# First the dimension of the river section is stated in whole positive numbers (integers)
# River has the following dimensions:
#   width = 15 centimeter until 200 centimeter
#   depth = 0 centimeter until 45 centimeter below watersurface
gridx = rng.randint(15, 200, size=1000)
gridy = rng.randint(0, 45, size=1000)
# Randomized integer numbers (cm) are now divided by 100 (m) and  depth is set to negative
# to match the original measured xy-dataset (m)
gridx = gridx/100
gridy = -gridy/100

# do kriging on unstructured grid
# -------------------------------
OK = gs.krige.Ordinary(model, [x, y], z, exact=True)
OK.unstructured([gridx,gridy])

# fast testplot for kriged data "river may be upside down"
# see  other code above

# get z-value from OK
krigz = OK.field.transpose()

# plot heatmap with plotly
# ----------------------------
# see  other code above