Variography#

import geolime as geo
from pyproj import CRS
import numpy as np
import pyvista as pv

pv.set_jupyter_backend('panel')


geo.Project().set_crs(CRS("EPSG:20350"))
/tmp/ipykernel_3153/3544329244.py:6: PyVistaDeprecationWarning: `panel` backend is deprecated and is planned for future removal.
  pv.set_jupyter_backend('panel')
dh = geo.read_file("../data/dh_pop_classif.geo")
dh.user_properties()
['X_COLLAR',
 'Y_COLLAR',
 'Z_COLLAR',
 'X_M',
 'Y_M',
 'Z_M',
 'X_B',
 'Y_B',
 'Z_B',
 'X_E',
 'Y_E',
 'Z_E',
 'Fe_pct',
 'Al2O3',
 'SiO2_pct',
 'K2O_pct',
 'CaO_pct',
 'MgO_pct',
 'TiO2_pct',
 'P_pct',
 'S_pct',
 'Mn_pct',
 'Fe_ox_ai',
 'hem_over_goe',
 'kaolin_abundance',
 'kaolin_composition',
 'wmAlsmai',
 'wmAlsmci',
 'carbai3pfit',
 'carbci3pfit',
 'Sample_ID',
 'Fe',
 'Fe2o3',
 'P',
 'S',
 'SiO2',
 'MnO',
 'Mn',
 'CaO',
 'K2O',
 'MgO',
 'Na2O',
 'TiO2',
 'LOI_100',
 'Depth',
 'ellipsoidal_distance',
 'OreZone',
 'domain_code',
 'domain',
 'manual_classif_code',
 'manual_classif',
 'sk_classif_code',
 'sk_classif']
dh.regions()
[]

Data Selection#

domain_solid = geo.datasets.load("rocklea_dome/domain_mesh.dxf")
domain_solid.name
'OreZone'
domain_solid.contains(dh)
domain_solid_pv = domain_solid.to_pyvista()
dh_pv = dh.to_pyvista('in_OreZone')
plotter = pv.Plotter()

plotter.add_mesh(domain_solid_pv, style="wireframe")
plotter.add_mesh(dh_pv.tube(radius=20))

plotter.set_scale(zscale=10)
plotter.show()
geo.histogram_plot(data=[{"object":dh, "property":"Fe_pct", "region":"in_OreZone"}])
dh.set_region_condition("HighGrade", "in_OreZone and (Fe_pct == Fe_pct)")

VarioMap#

lags, tol = geo.generate_lags(lag=200, plag=50, nlags=15)
geo.vario_map(
    geo_object=dh,
    attribute="Fe_pct",
    region="HighGrade",
    lags=lags,
    tol=tol,
    n_az=15,
    atol=35
)
geo.vario_contour(
    geo_object=dh,
    attribute="Fe_pct",
    region="HighGrade",
    lags=lags,
    tol=tol,
    n_az=15,
    atol=35
)
../_images/fed0f72945d5227ca591f7dbefac5b2162ff1486bb444cfd16110eb24a6e1fab.png

VarioMap suggets major and minor preferential direction are N00 and N90.

major_direction = 0.
minor_direction = 90.

Fitting a Covariance Model#

Model the Nugget Effect#

Experimental Downhole Variogram#

lags_bh, tol_bh = geo.generate_lags(lag=1., plag=50, nlags=10.)
tol_bh
0.5
vario_exp_bh = geo.variogram(
    object=dh,
    attribute="Fe_pct",
    region="HighGrade",
    geographic_azimuth=0,
    dip=90,
    pitch=90,
    lags=lags_bh,
    tol=tol_bh,
    atol=10
)

geo.plot_semivariogram(
    variograms=[vario_exp_bh],
    display_npairs=True
)
vario_exp_bh
rank npairs lag avgdist vario indices
0 1 1601 1.0 1.0 38.453467 [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6...
1 2 1530 2.0 2.0 60.226471 [[0, 2], [1, 3], [2, 4], [3, 5], [4, 6], [5, 7...
2 3 1460 3.0 3.0 71.746233 [[0, 3], [1, 4], [2, 5], [3, 6], [4, 7], [5, 8...
3 4 1394 4.0 4.0 80.208393 [[0, 4], [1, 5], [2, 6], [3, 7], [4, 8], [5, 9...
4 5 1324 5.0 5.0 83.884819 [[0, 5], [1, 6], [2, 7], [3, 8], [4, 9], [5, 1...
5 6 1256 6.0 6.0 83.457803 [[0, 6], [1, 7], [2, 8], [3, 9], [4, 10], [5, ...
6 7 1192 7.0 7.0 84.465604 [[0, 7], [1, 8], [2, 9], [3, 10], [4, 11], [5,...
7 8 1124 8.0 8.0 82.919039 [[0, 8], [1, 9], [2, 10], [3, 11], [4, 12], [5...
8 9 1057 9.0 9.0 83.146168 [[0, 9], [1, 10], [2, 11], [3, 12], [4, 13], [...

Automatic Fitting#

nugget_model = geo.Nugget() + geo.Spherical()
geo.model_fit(
    variograms=[vario_exp_bh], 
    cov=nugget_model
)
geo.plot_semivariogram(
    variograms=[vario_exp_bh],
    model=nugget_model,
    model_angles=[{"azi":0, "dip":90, "pitch":90}],
    display_npairs=True
)
print(nugget_model)
Model with 2 components 
 
Component 1 : 
Sill : 16.178499493157535 
Covariance type : Nugget 

Component 2 : 
Sill : 65.8248169561377 
Covariance type : Spherical 
Scales : (7.577161623484672, 2.25, 4.286255652839992) 
Angles : (1e-10, 1e-10, 0.0) 

Total sill = 82.00331644929523 
nugget_value = nugget_model.cov_elem_list[0].sill

Model the Sphericals#

Experimental Variograms along the XY plane#

lags, tol = geo.generate_lags(lag=50, plag=50, nlags=10)
vario_exp_major = geo.variogram(
    object=dh,
    attribute="Fe_pct",
    region="HighGrade",
    geographic_azimuth=major_direction,
    dip=0,
    pitch=0,
    lags=lags,
    tol=tol,
    atol=40
)
vario_exp_minor = geo.variogram(
    object=dh,
    attribute="Fe_pct",
    region="HighGrade",
    geographic_azimuth=minor_direction,
    dip=0,
    pitch=0,
    lags=lags,
    tol=tol,
    atol=30
)

geo.plot_semivariogram(
    variograms=[vario_exp_major, vario_exp_minor],
    display_npairs=True
)

Automatic Fitting#

model =  geo.Nugget() + geo.Spherical() + geo.Spherical()
short_range_spherical_ratio = 0.7
long_range_spherical_ratio = 0.3
geo.model_fit(
    variograms=[
        vario_exp_major, 
        vario_exp_minor, 
        vario_exp_bh
    ], 
    cov=model,
    constraints=[
        {"sill_fixed":nugget_value},
        {
            "sill_max": short_range_spherical_ratio * (vario_exp_major.attrs["variable_var"] - nugget_value),
            "angle_fixed_0": major_direction,
            "angle_fixed_1": 0,
            "angle_fixed_2": 0
        },
        {
            "sill_max": long_range_spherical_ratio * (vario_exp_major.attrs["variable_var"] - nugget_value),
            "angle_fixed_0": major_direction,
            "angle_fixed_1": 0,
            "angle_fixed_2": 0,
        }
    ]
)
geo.plot_semivariogram(
    variograms=[
        vario_exp_major, 
        vario_exp_minor, 
        vario_exp_bh
    ],
    model=model,
    model_angles=[{"azi":major_direction, "dip":0, "pitch":0}, {"azi":minor_direction, "dip":0, "pitch":0}, {"azi":0, "dip":90, "pitch":90}],
    display_npairs=True
)
print(model)
Model with 3 components 
 
Component 1 : 
Sill : 16.178499493157535 
Covariance type : Nugget 

Component 2 : 
Sill : 53.16450791195954 
Covariance type : Spherical 
Scales : (199.95667759111325, 66.21593638344216, 3.8237237637287356) 
Angles : (0.0, 0.0, 0.0) 

Component 3 : 
Sill : 22.784789044789914 
Covariance type : Spherical 
Scales : (91.91627353545513, 85.32983957687435, 14.111597684097493) 
Angles : (0.0, 0.0, 0.0) 

Total sill = 92.127796449907