RMSP Highlights


Getting Started with Kriging in RMSP

Ryan Barnett,

RMSP includes many example data sets for testing algorithms. Here, a small underground copper data set is used to introduce kriging in RMSP. Start by importing RMSP and loading the example data set and summarizing it with some basic statistics. Here, we only keep the drill hole ID, composite centroid, and copper assay.

import rmsp
comps = rmsp.load_example_data('copper_composite')[["DHID", "x", "y", "z", "CU"]]
comps.describe()
x y z CU
count1310 1310 1310 499
mean 872.62312259 1148.67 2.66148
std 140.019 46.1881 71.9665 2.31973
min 615.52612125.6 932.47 0
25% 743.53812226.3 1099.75 0.795872
50% 898.89612258.8 1155.41 2.17811
75% 982.17812291.1 1202.22 3.91401
max 1109.04 12409.1 1310.79 14.1608

Copper grades (CU) are estimated here. A dictionary is used for storing the limits (clim) and log scaling (cmap_log) that is used for coloring copper in section plotting. Similarly, a dictionary stores spatial limits and figure size that is used for sections.

cu_pars = dict(clim=(1, 10), cmap_log=True)
plot_pars = dict(figsize=(12, 12), ylim=(12175, 12375), xlim=(550, 1250))
comps.sectionplot('CU', s=20, **plot_pars, **cu_pars)
Infinite section of example copper grades.

Load a triangulated solid from an STL file, which bounds a copper estimation domain that will be applied. The resulting solid object is rmsp.Solid type, providing functionality for working with triangulations.

solid = rmsp.load_example_data('copper_low_shell')

Define a coordinate for slicing through the solid, composites and other objects in this notebook. Define tolerances that control the viewing corridor. Also define some properties specific to composites and solids.

coordinate = 1100
solid_pars = dict(coordinate=coordinate, tolerances=5, face_c='C0',
                  edge_c=None, face_label='Solid')
comps_pars = dict(coordinate=coordinate, tolerances=25, s=20,
                  label='Composites')

Compare the composites and solids in sections with pre-defined properties. Apply the plotting in a function so that we can call this again without repeating code (after modifying comps).

def plot_comps_solid():
    fig, ax, cax = comps.sectionplot('CU', **cu_pars, **comps_pars, **plot_pars)
    solid.sectionplot_draw(ax, **solid_pars)
    ax.legend()

plot_comps_solid()
Infinite section of example copper grades and a copper grade shell.

Reduce the data to records that will be used for estimation. Flag comps records inside solid using rmsp.Solid.contains, before reducing comps to those records using the standard PointData.loc[] slicing. Further reduce comps by removing null records.

insolid = solid.contains(comps)
print(f'there are {insolid.sum()} of {len(comps)} records inside the solid')
comps = comps.loc[insolid].reset_index(drop=True)
comps = comps.loc[~comps['CU'].isna()].reset_index(drop=True)
print(f'there are {len(comps)} non-null records that remain for estimation')

# Outputs:
"""
there are 598 of 1310 records inside the solid
there are 239 non-null records that remain for estimation
"""

Plot the updated comps against the wireframe using the previously defined function.

plot_comps_solid()
Infinite section of example copper grades clipped to grade shell.

The rmsp.GridData is the standard object for storing and working with grid data. Dimensions of the grid are \(uvz\) rather than \(xyz\), since an azimuth rotation may be used in the definition. A large variety of methods exist for the manual and automated generation of GridData, including automated population of an rmsp.Solid as used here.

grid = rmsp.GridData.from_solid(solid, usize=15.0, vsize=7.5, zsize=7.5,
                                azimuth_deg=-15.0, srch_radius=50.0,
                                origin=(600.0, 12000.0, 930.0))

Define grid plotting parameters before comparing against the wireframe. Note that rmsp.GridData is sparsely populated - its storage and the memory requirements of all RMSP routines in which it is applied are a function of populated blocks.

grid_pars = dict(block_outline=True, coordinate=coordinate, label='Grid')

fig, ax, cax = grid.sectionplot('.6', **grid_pars, **plot_pars)
solid.sectionplot_draw(ax, **solid_pars)
ax.legend()
Grid generated inside copper grade shell.

Define variable lag spacing with rmsp.Lags; ten 3m lags with 3m tolerances, followed by five 10m lags with 10m tolerances. Then define an omni-directional experimental variogram search (ExpVarioSearch), where the lags described above are calculated in equal spaced directions. The associated experimental variograms are then calculated with rmsp.ExpVario().calculate.

lags = [rmsp.Lags(3.0, 3.0, 10), rmsp.Lags(10.0, 10, 5)]
expsearches = rmsp.ExpVarioSearch(0.0, 0.0, lags, azmtol=90, incltol=90)
expvario = rmsp.ExpVario().calculate(comps, 'CU', expsearches)

Fit a variogram model rmsp.VarioModel to the experimental variogram before plotting the results. Also output the variogram parameters in a table.

variomodel = rmsp.VarioModel.fit_experimental(expvario, nugget=0.1, num_struct=2,
                                              ranges12_bounds=1, ranges13_bounds=1)
fig, ax = expvario.plot()
variomodel.plot_draw(ax)
variomodel.to_table()
Omnidirectional copper variogram model.
Nugget Structure 1 Structure 2
Contribution0.100 0.309 0.591
Model Shape spherical spherical
Angle 1 0.0 0.0
Angle 2 0.0 0.0
Angle 3 0.0 0.0
Range 1 154.0 5.0
Range 2 154.0 5.0
Range 3 154.0 5.0

Setup an estimation search ellipse (rmsp.Search) with principal angles of 86°, 8° dip and 0°, ranges of 200x50x200m and typically used min/max constraints.

search = rmsp.Search([86.0, 8.0, 0.0], [200.0, 50.0, 200.0],
                     min_comps=2, max_comps=8, max_comps_per_dhid=4,
                     max_comps_per_oct=5)

Setup a kriging estimator (rmsp.KrigeEstimator), which in this case performs block ordinary kriging using 2x3x2 discretization, before estimating across the grid.

estimator = rmsp.KrigeEstimator().set_params(search, variomodel, 'ok', num_discret=[2, 3, 2])
ests = estimator.estimate(grid, comps, 'CU', num_threads=5)

Plot sections of the estimates with the data and wireframes.

fig, ax, cax = ests.sectionplot('estimate', **cu_pars, **grid_pars, **plot_pars)
comps.sectionplot_draw(ax, 'CU', edgecolors='k', **cu_pars, **comps_pars)
solid.sectionplot_draw(ax, **solid_pars)
ax.legend()
Section of kriged copper estimates.

Interested in using RMSP for your resource modeling challenges?

  contact@resmodsol.com