anisotropy angle definitions, pykrige vs gstools

[collord]

I'm new and trying to wrap my head around the angle definitions for anisotropy. There seems to be a difference between how pykrige and gstools handle things. Here is a simple ellipsoid visualization that looks like dip (rotation about X) is doubling the angle in the non-gstools variogram, but correct in gstools version. Or tell me where I've gone wrong! Code below.

With gstools:

Without gstools:

from pykrige.ok3d import OrdinaryKriging3D

import numpy as np
import gstools as gs
import pyvista as pv
gridx = np.arange(0.0, 1.0, 0.01)
gridy = np.arange(0.0, 1.0, 0.01)
gridz = np.arange(0.0, 1.0, 0.01)
data = np.array(
    [
        [0.1, 0.1, 0.1, 0.5],
        [0.9, 0.1, 0.1, 0.5],
        [0.1, 0.1, 0.9, 0.5],
        [0.1, 0.9, 0.1, 0.5],
        [0.9, 0.9, 0.1, 0.5],
        [0.1, 0.9, 0.9, 0.5],
        [0.9, 0.1, 0.9, 0.5],
        [0.9, 0.9, 0.9, 0.5],
        [0.5, 0.5, 0.5, 1.0],
    ]
)
dips = [float(i)*(np.pi/8.0) for i in range(5)]
dipdirs = [float(i)*(np.pi/4.0) for i in range(2)]
use_gs = False
p = pv.Plotter(window_size=(1400, 2000), shape=(len(dips), len(dipdirs)) )
for i,dip in enumerate(dips):
    for j,dipdir in enumerate(dipdirs):
        if use_gs:
            ok3d = OrdinaryKriging3D(
                data[:, 0], data[:, 1], data[:, 2], data[:, 3],
                gs.Gaussian(dim=3, len_scale=[0.05,0.4,0.2], angles=[dip,0,dipdir], var=0.1, nugget=0.0) 
            )            
        else:
            ok3d = OrdinaryKriging3D(
                data[:, 0], data[:, 1], data[:, 2], data[:, 3], 
                variogram_model="gaussian", variogram_parameters ={'psill': .1, 'range': .05, 'nugget': .0},
                anisotropy_angle_x=dipdir*(360.0/(2.0*np.pi)), anisotropy_angle_y=0.0, anisotropy_angle_z=dip*(360.0/2.0*np.pi),
                anisotropy_scaling_y=0.1,anisotropy_scaling_z=.2
            )
        k3d1, ss3d = ok3d.execute("grid", gridx, gridy, gridz)
        grid = pv.UniformGrid()
        grid.dimensions = np.array(k3d1.shape) + 1
        grid.origin = (0, 0, 0)  
        grid.spacing = (.01, .01, .01)  
        grid.cell_arrays["values"] = k3d1.flatten(order="F")  
        p.subplot(i,j)
        p.add_volume(grid,opacity = [0, 0, 0, 0, 0, 1, 1], shade=True)
        p.add_text('dipdir%.2f dip%.2f '%( dipdir,dip), font_size=12, shadow=True)
        p.show_grid() 
p.show(cpos=[(0,0,3),(0.5,0.5,0.5),(0,1,0)])

Thanks

[MuellerSeb]

Here I am meeting my nemesis again. But I think there is nothing wrong. The rotation angle vector in GSTools is sorted by rotation planes called Tait–Bryan angles or yaw, pitch and roll:

1. xy-plane angle (corresponds to anisotropy_angle_z in PyKrige) -> yaw
2. xz-plane angle (corresponds to anisotropy_angle_y in PyKrige) -> pitch
3. yz-plane angle (corresponds to anisotropy_angle_x in PyKrige) -> roll

So, in your case this should be fixed by swapping anisotropy_angle_x and anisotropy_angle_z.

Hope that helps! Seb

[MuellerSeb]

@collord did that solve your problem?

[collord]

@MuellerSeb Yes, I think that makes sense and results look right. If the projects are interrelated now, might make sense to standardize on the angle definitions? Thank you again for the hint.

[MuellerSeb]

@collord we will do so in the future, since the GSTools CovModel will be the standard way to provide the model.