NumPy Error when Ordinary Kriging

[banesullivan]

I'm trying to teach myself how to use this library and thought I'd try out oridinary kriging but I keep running into an issue when I try to krige a random sample of points from a known dataset:

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

# Generate some test data
original = examples.load_channels().ctp()
n = original.n_points // 5
print("{:.1f}%".format(n / original.n_points * 100))
idx = np.random.randint(0, original.n_points, size=n)
cond_pos = original.points[idx]
cond_val = original['facies'][idx]

# conditions
model = Gaussian(dim=3,)
krig = krige.Ordinary(model, cond_pos=cond_pos, cond_val=cond_val)

# create the kriging grid using the original dataset
x = np.unique(original.points[:, 0])
y = np.unique(original.points[:, 1])
z = np.unique(original.points[:, 2])

# Krige it - ERROR OCCURS HERE
krig((x,y,z), mesh_type="structured")

# Grab a PyVista mesh
krigged = krig.to_pyvista()
krigged.plot()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-d204afd9cf26> in <module>
     22 
     23 # Krige it
---> 24 krig((x,y,z), mesh_type="structured")
     25 
     26 # Grab a PyVista mesh

~/Documents/OpenGeoVis/External/GSTools/gstools/krige/ordinary.py in __call__(self, pos, mesh_type)
    111             (np.full_like(self.cond_val, self.model.sill), [1])
    112         )
--> 113         self.mean = np.einsum("i,ij,j", cond, krig_mat, mean_est)
    114 
    115         # reshape field if we got an unstructured mesh

~/anaconda3/envs/dev/lib/python3.7/site-packages/numpy/core/einsumfunc.py in einsum(*operands, **kwargs)
   1344     # If no optimization, run pure einsum
   1345     if optimize_arg is False:
-> 1346         return c_einsum(*operands, **kwargs)
   1347 
   1348     valid_einsum_kwargs = ['out', 'dtype', 'order', 'casting']

ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (1272621,)->(1272621,newaxis) (4,4)->(4,4) (1272621,)->(1272621) 

[MuellerSeb]

Hey there,
you made a small mistake in your script. The cond_pos needs to be transposed in your case, since the first dimension of this array should refer to the x-, y- and z-axis.

So all you have to change is:

cond_pos = original.points[idx].T

Cheers,
Sebastian

[LSchueler]

Added a readable exception message for such an input with commit efe78a7.

[banesullivan]

Thanks for the feedback!

Is using the first axis as the x,y,z axes consistent across the API? I didn't have an error when experimenting with Simple Kriging (however my results weren't great so maybe that's why)

[banesullivan]

Ahh actually yes, I'm seeing now that this is consistent across the examples.

[banesullivan]

Also, that dataset above was taking forever, so I'm trying to use a much smaller dataset now.

I'm curious if you all could help me create a few examples of 3D kriging. Perhaps we start with this example from PyVista where we try kriging some 3D point data onto a "terrain" surface in 3D (in pyvista, we simply use a Gaussian kernel based interpolation)? This is what I have so far, but it's not really working:

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

surface = examples.download_saddle_surface()
samples = examples.download_sparse_points()

p = pv.Plotter()
p.add_mesh(samples, point_size=30.0, render_points_as_spheres=True)
p.add_mesh(surface)
p.show()

cond_pos = samples.points.T
cond_val = samples['val']

# conditions
model = Gaussian(dim=3,)
krig = krige.Ordinary(model, cond_pos=cond_pos, cond_val=cond_val)

# create the kriging grid using the original dataset
x = surface.points[:, 0]
y = surface.points[:, 1]
z = surface.points[:, 2]

# Krige it
krig((x,y,z), mesh_type="structured")

# Grab a PyVista mesh
krigged = krig.to_pyvista()
krigged.plot()

... ran for a while then killed my kernel. My system info is as follows:

>>> pv.Report(["scipy", "hankel", "emcee", ])
--------------------------------------------------------------------------------
  Date: Wed Sep 25 17:02:15 2019 MDT

            Darwin : OS
                12 : CPU(s)
            x86_64 : Machine
             64bit : Architecture
           32.0 GB : RAM
           Jupyter : Environment

  Python 3.7.3 | packaged by conda-forge | (default, Jul  1 2019, 14:38:56)
  [Clang 4.0.1 (tags/RELEASE_401/final)]

             1.3.0 : scipy
             0.3.9 : hankel
            3.0rc2 : emcee
            0.22.2 : pyvista
             8.2.0 : vtk
            1.16.3 : numpy
             2.5.0 : imageio
             1.4.3 : appdirs
             0.4.3 : scooby
             3.1.0 : matplotlib
             5.9.2 : PyQt5
             7.6.1 : IPython
             7.4.2 : ipywidgets
             1.0.0 : colorcet
               2.0 : cmocean
             0.6.2 : panel

  Intel(R) Math Kernel Library Version 2018.0.3 Product Build 20180406 for
  Intel(R) 64 architecture applications
--------------------------------------------------------------------------------

[banesullivan]

Ahh... turns out I needed an "unstructured" mesh type. Doing so and using some array ordering majic, I can get very similar results as the example in PyVista!

cond_pos = samples.points.T
cond_val = samples['val']

# conditions
model = Gaussian(dim=3, len_scale=12)
krig = krige.Ordinary(model, cond_pos=cond_pos, cond_val=cond_val)

# create the kriging grid using the original dataset
x = surface.points[:, 0]
y = surface.points[:, 1]
z = surface.points[:, 2]

# Krige it
krig((x,y,z), mesh_type="unstructured")

# Grab a PyVista mesh
surface["krige"] = krig.field
surface.plot(scalars="krige")

[banesullivan]

So this leaves me with one final question, are "structured" grids, strictly rectilinear? I.e the surface I use above is not rectilinear and thus "unstructured"? Solving it this way is totally fine as PyVista can manage all of the geometry and all GSTools needs are the node locations...

[LSchueler]

Sure we can help! - But obviously, you're figuring out by yourself ;-)

You're now the second one to stumble upon the nomenclature of "structured" / "unstructured".
You are right, "structured" can only be used for rectilinear grids, as in it can be spanned by X individual arrays, where X is the spatial dimension.
For everything else you have to provide a position tuple for every field value.
The "structured" grid saves some memory and maybe more importantly, it can make more abstract applications way easier to handle, e.g. if something is aligned with one of the axes.

What exactly did you think the "structured" mesh would be?

[LSchueler]

By the way, in case you really have trouble with long computation times, have you tried out compiling GSTools with OpenMP support? - It makes quite a difference:

pip install --global-option="--openmp" https://github.com/GeoStat-Framework/GSTools/archive/master.zip