Is it possible with the current implementation to conduct Kriging for unstructured grid or a set of (x, y) points?

[flydream0428]

In the example/tutorial notebook, Kriging is always done for a structured regular-spaced grid. I am wondering if it is possible to only conduct kriging for a list of (x, y) points? Thanks!

[PauloCarvalhoRJ]

Hello,

Depending on your application, you can model the variogram and perform kriging in depositional regular space (UVW) instead of Cartesian space (XYZ). The figure below illustrates an irregular grid representing some folded geology (A) and the corresponding regular grid in UVW space (B):

It's easy to follow but there's a caveat: you have to unfold the primary data (e.g. well or borehole data), as illustrated in the figure below:

So, it would be necessary for GeostatsPy to include an unfolding routine as described here: https://github.com/PauloCarvalhoRJ/gammaray/blob/master/docs/GammaRayManual.docx (Section 15.3) and implemented in C++ here: https://github.com/PauloCarvalhoRJ/gammaray/blob/master/domain/geogrid.cpp (look for function GeoGrid::unfold()). Of course, there are other unfolding algorithms out there. It's likely there is one already implemented in Python.

regards,

Paulo

[PauloCarvalhoRJ]

The figure below illustrates how to compute 2D stratigraphic coordinates of a sample point P with respect to an irregular cell of a 2D unstructured surface:

source: https://math.stackexchange.com/questions/13404/mapping-irregular-quadrilateral-to-a-rectangle

This computation is used to compute texture coordinates in computer graphics. So, I extended it to 3D to compute stratigraphic coordinates in irregular meshes (reference to the code posted above).

[flydream0428]

@PauloCarvalhoRJ Thanks a lot for your reply, Paulo. I may be wrong here, but isn't easier to just add a function to allow Kriging on a list of points based on their locations rather than a structured grid? When we supply a structured grid to GSLIB, my understanding is it will calculate the coordinates internally anyway. The function can return the results as a list rather than a matrix.

[PauloCarvalhoRJ]

One of the reasons behind using a regular grid is efficiency. But if you wish to krige over a point set, you can study the kb2d.py's kb2d() function code to create a new function intended for point sets. Here's the loop over the grid blocks: https://github.com/GeostatsGuy/GeostatsPy/blob/a256e1b271b33a413ede8eb3761b92ac7aa4059b/working/kb2d.py#L109-L239. You'd need to modify that loop to read the explicit XYZ coordinates coming from the point set instead of computing them from the grid specs. You also have to modify the part dealing with the block variance, afterall you'd be kriging on points.