Ordinary Kriging: help please?

[banesullivan]

So I'm trying to get up and running with this software, and I'm finding it rather difficult to produce a dirty ordinary kriging result. I learned geostats/kriging using SGeMS and I'm realizing much of what I learned is difficult to transfer here (so far)... I'm curious if someone could help me get some quick/dirty results with a dataset that I've used in SGeMS before.

The Example

I believe my Gaussian model is probably problematic

import numpy as np
from gstools import Exponential, Gaussian, krige
import pyvista as pv

# Load the data
probes = pv.read("probes.vtp")
model_space = pv.read("model_space.vtr")

model = Gaussian(dim=3, var=10, len_scale=500)
cond_pos = probes.points.T # gstools convention
# conditional values - what i want to interpolate
field = probes['temperature (C)']
# Create the kriging model
krig = krige.Ordinary(model, cond_pos=cond_pos, cond_val=field)

krige_field, krige_var = krig((model_space.points.T))

model_space["krige_field"] = krige_field
model_space["krige_var"] = krige_var

model_space.contour().plot(scalars="krige_field" )

I think I am running into trouble with the len_scale option in the Gaussian model. GSTools enforces a limit of 1000 on that but for this dataset, I may need a higher len_scale?

If you all were to use this given dataset to produce an example showcasing what GSTools can do, what would that be? This is all in an effort to help me build out how PyVista might interface with GSTools

Data

I have sparse, subsurface temperature probe data and want to krige a volume of these data.

data.zip

"Truth"

I'm working towards results that look more like this (prodcued in SGeMS):

[LSchueler]

Without having had a more in depth look at your problem yet and if the length scale bound is the problem, you can do following for example:

model = Gaussian(dim=3, var=10)
model.set_arg_bounds(len_scale=[0, 10000])
model.len_scale = 2000

[LSchueler]

If you pull the newest version of GSTools, the boundaries are now much generous ;-)

[banesullivan]

Okay, awesome! I pulled the latest changes and increasing the length scale is helping me approach a better solution.

Do you all have any examples, tips, or tricks for finding a good Gaussian model?

[javoha]

Hello,

let's see if I can make a suggestion (please correct me if I got anything wrong):

The len_scale of GSTools corresponds to the range describing the spatial correlation of data when performing variogram analysis. There is basically two ways to obtain a range for your model.

1. Expert knowledge: You have some idea of the length scale of spatial correlation in your data set that you use as a range. In this case using the same range (len_scale) that you used for the SGeMS model should give you the same result. Just be aware that some software packages use the "effective range" (for the Gaussian model this is sqrt(3)*range I think).

2. You perform a variogram analysis (basically a data driven approach to obtain a range). If I am correct GSTools offers this: https://geostat-framework.readthedocs.io/projects/gstools/en/latest/tutorial_03_vario.html

Apart from that you might have anisotropies in your model (in your example it looks like you have strong continuity in x and y direction (horizontal) and variation in z direction corresponding to the direction of sedimentation. Also in literature (e.g. Webster and Oliver (2008)) the Gaussian model is deprecated for approaching the origing with zero gradient. You might want to consider a different variogram model.

Maybe this is helpful and please correct me if I got anything wrong.

[LSchueler]

Hi @banesullivan,

did the reply from @javoha (thanks for that!) help? - Could you reproduce the "truth" from SGeMS?

[banesullivan]

Thanks for the tips @javoha!

I haven't had the time to jump back into this... I will need to next week-ish and I will report back!

You might want to consider a different variogram model.

@javoha, Are there any other options for models when kriging with GSTools?

@LSchueler, feel free to close this issue as it's not really an issue/bug with the software but a support topic

[MuellerSeb]

@banesullivan : you can use all provided models for kriging: provided-covariance-models
And of course self defined ones. ;-)

[MuellerSeb]

Closing now. But thanks for the input ;-)

[banesullivan]

So I'm trying this again but with another dataset that is a bit cleaner. Data attached and this data has a permissive license if you want to use it in an example in your docs - just let me know.

I've compiled with OpenMP support per #31 (comment) but as soon as I start running the kriging (either on a structured or unstructured mesh) the memory consumption blows up with over 100 GB used on the SWAP so I kill the task. This amount of memory usage seems ridiculous to me...

@MuellerSeb, @LSchueler, or @javoha: Can you all help in some capacity to make this example run?

Or can you put some 3D examples in the docs besides just a random field? I can't get anything to work in 3D with this library yet.

Data files: Archive.zip

import numpy as np
import matplotlib.pyplot as plt

import pyvista as pv

from gstools import (Exponential, Gaussian, krige, 
                     vario_estimate_unstructured, Stable)

grid = pv.read("grid.vtr")
assay = pv.read("assay.vtp")

cond_pos = assay.points.T # gstools convention
# conditional values - what i want to interpolate
field = assay["CU_pct"]


model = Exponential(dim=3, var=2, len_scale=10)
# estimate the variogram of the field with 40 bins
bins = np.arange(40)
bin_center, gamma = vario_estimate_unstructured(cond_pos, field, bins)
# fit the variogram with a stable model. (no nugget fitted)
fit_model = Stable(dim=3)
fit_model.fit_variogram(bin_center, gamma, nugget=False)
# output
ax = fit_model.plot(x_max=40)
ax.plot(bin_center, gamma)
print(fit_model)

Not too shabby of a variogram, so I should be able to at least produce some results with this....

# Create the kriging model
krig = krige.Ordinary(model, 
                      cond_pos=cond_pos, 
                      cond_val=field)

And this is where the memory usage blows up

krige_field, krige_var = krig((grid.x, grid.y, grid.z), mesh_type="structured")

kgrid["CU_pct"] = krige_field
kgrid["krige_var"] = krige_var

[banesullivan]

Or as unstructured:

krige_field, krige_var = krig(grid.points.T, mesh_type="unstructured")

kgrid["CU_pct"] = krige_field
kgrid["krige_var"] = krige_var

[banesullivan]

And...

model = Exponential(dim=3, var=2, len_scale=10)
# estimate the variogram of the field with 40 bins
bins = np.arange(100)
bin_center, gamma = vario_estimate_unstructured(cond_pos, field, bins)
# fit the variogram with a stable model. (no nugget fitted)
fit_model = Stable(dim=3)
fit_model.fit_variogram(bin_center, gamma, nugget=False)
# output
ax = fit_model.plot(x_max=100)
ax.plot(bin_center, gamma)
print(fit_model)