Issues with Kriging values

[donhalmina]

I am performing an ordinary kriging and the results I get after performing the kriging are totally wrong. For simplicity, my grid points are also my data points. The working code is:

"Import libraries and modules"
import numpy as np
import gstools as gs
import tkinter as tk
from tkinter import filedialog
import pandas as pd
from scipy.spatial.distance import pdist, squareform

"Import data"
root = tk.Tk()
root.withdraw()
filepath = filedialog.askopenfilenames(parent=root, title='Choose a file')
data = pd.DataFrame()
for file in filepath:
    df = pd.read_excel(file)
    data = data.append(df)
data = data.iloc[:, :3]
data = data[~np.isnan(data).any(axis=1)]
data = np.array(data, dtype=np.float)

"Calculate semivariance"
sill = np.var(data[:, 2])
pwdist = squareform(pdist(data[:, :2]))
max_lag = np.int(np.ceil(np.max(pwdist)))
npoint = np.int(np.floor(max_lag/40)) 
firstlag = max_lag / 100
lags = np.linspace(firstlag, max_lag, npoint)
bin_center, gamma = gs.vario_estimate_unstructured((data[:,0], data[:,1]), data[:,2], lags)

"Fit theoretical variogram to data"
fit_model = gs.Stable(dim=2)
fit_model.set_arg_bounds(var=[0, sill])
fit_model.set_arg_bounds(len_scale=[0, 1000000])
fit_model.fit_variogram(bin_center, gamma, nugget=False)
ax = fit_model.plot(x_max=max_lag)
ax.plot(bin_center, gamma)
print(fit_model)

"Perform kriging"
OK = gs.krige.Ordinary(fit_model, [data[:, 0], data[:, 1]], data[:, 2])
OK.unstructured([data[:, 0], data[:, 1]])
ax = OK.plot()
ax.set_aspect("equal")

The semivariogram looks good (as attached below).

The value of my kriging should range between 0 and 400 as data[:, 2] ranges between this value. However, i am getting values between -3000 and 8000.

Find attached the data to reproduce the problem

CNDA-ESP_SIMRAD_MDS_DGPS_2000_2001.xlsx

[MuellerSeb]

What helped in my case, was to change the estimator to cressie:

bin_center, gamma = gs.vario_estimate_unstructured(
    (data[:, 0], data[:, 1]), data[:, 2], lags, estimator="cressie"
)

But the only difference is a different estimated len-scale. So I think the problem is in the kriging routines. An educated guess is, that the resulting kriging matrix is singular and produces an invalid inverse matrix. Maybe this will be resolved by this suggestion: #79 (comment)

I will keep track of that!

[donhalmina]

I changed the estimator to Cressie, to try to reproduce your results but still get erreoneous result. I am using gstools version 1.2.1 and don't know if it is the same version you are using as I am expecting to obtain the same results you obtianed.

I have my Geostat tool which i programmed and use and for determining the weight parameter, I use np.linalg.lstsq as it is more robust and stable than np.linalg.inv.

Also, how can I obtain the values of the kriging estimate or export it (.csv, .txt, etc)?

[MuellerSeb]

We will use the pseudo-inverse in the future, which should hold equal results as using least-square.

Exporting can be done to a "vtk" file as demonstrated here:
https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/00_misc/01_export.html#sphx-glr-examples-00-misc-01-export-py
and here:
https://geostat-framework.readthedocs.io/projects/gstools/en/stable/index.html#vtk-pyvista-export

[MuellerSeb]

This problem will be solved with #95

[MuellerSeb]

This should be solved with #97

Feel free to re-open this issue, if you have further questions.