RMSP Highlights

Getting Started with Conditional Simulation 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. Load composited copper data using the load_example_data function, returning a comps object that is rmsp.PointData type. Reduce the composites to records with measured copper.

import rmsp
comps = rmsp.load_example_data('copper_composite')[["DHID", "x", "y", "z", "CU"]]
comps = comps.loc[~comps['CU'].isna()].reset_index(drop=True)
DHID x y z CU
0DDH-0001820.76712272.91133.6 3.16213
2DDH-0001820.9 12276.41124.243.87915
4DDH-0001821.07 12279.91114.872.78315

Generate a grid for simulation using rmsp.GridData.from_scattered_points. Grid nodes are generated in this case at locations within 50m of composite data.

grid = rmsp.GridData.from_scattered_points(
    comps, usize=5.0, vsize=2.5, zsize=2.5, azimuth_deg=-15.0,
    srch_radius=50.0, origin=(540.0, 12050.0, 930.0))
print(f'the grid has {len(grid)} nodes')

# Outputs
the grid has 532553 nodes

Define a slicing coordinate that is used throughout this example, as well as plotting parameters for copper and each object within respective dictionaries.

coordinate = 1110
cu_pars = dict(clim=(1, 10), cmap_log=True)
comps_pars = dict(coordinate=coordinate, tolerances=5, s=20, edgecolors='k',
grid_pars = dict(block_outline=True, coordinate=coordinate, label='Grid')
plot_pars = dict(figsize=(12, 12), xlim=(600, 1100), ylim=(12150, 12425))

Compare the composites and generated grid in a horizontal section using the dictionaries defined above.

fig, ax, cax = comps.sectionplot('CU', **cu_pars, **comps_pars, **plot_pars)
grid.sectionplot_draw(ax, '.6', **grid_pars)
Section of grid generated around copper composites.

Simulation requires normal score conditioning data and its associated variogram model. Begin by despiking the data with a spatial moving average. Otherwise random breaking of spikes are applied internally by the normal score transform.

despike = rmsp.DespikeSpatial().fit(comps['CU'])
print(f'there are {despike.spike_indices.sum()} records with spikes')
comps['CU'] = despike.transform(comps, 'CU')

# Outputs
there are 21 records with spikes

Decluster the data using nearest neighbour weighting (NN). These weights are used for applying the normal score transform and checking simulation results. Alternative declustering methods are available in RMSP, such cell, ordinary kriging and trend based.

comps['wt'] = rmsp.DeclusterNN([0.0] * 3, [1000.0] * 3).decluster(grid, comps)

Compare statistics of the clustered and representative distributions.

fig, axes = plt.subplots(1, 2, figsize=(8, 4), sharey=True)
comps.histplot('CU', bin_type='fraction', ax=axes[0])
comps.histplot('CU', 'wt', bin_type='fraction', ax=axes[1])
Comparison of copper statistics between clustered and declustered data.

Apply the normal score transform in preparation for Gaussian simulation. The nscore object of rmsp.NSTransformer type stores the transform table that is required for back-transforming realizations.

nscore = rmsp.NSTransformer()
comps['NS_CU'] = nscore.fit_transform(comps['CU'], comps['wt'])
fig, ax = comps.histplot('NS_CU', 'wt', histlim=[-3, 3], xlim=[-3, 3])
Normal score transformed distribution.

Calculate an experimental normal scores variogram with an omni-directional search and variable lag spacing.

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, 'NS_CU', expsearches)

Fit a variogram model rmsp.VarioModel to the experimental variogram before plotting the results.

variomodel = rmsp.VarioModel.fit_experimental(expvario, nugget=0.1, num_struct=2,
                                          ranges12_bounds=1, ranges13_bounds=1)
fig, ax = expvario.plot()
Normal scored copper omnidirectional variogram.

Set the number of realizations to simulate and the search ellipse parameters rmsp.Search, before performing conditional simulation using rmsp.Simulator. The realizations are stored in an efficient binary format that is accessed via the returned simcache of rmsp.SimCache type.

num_reals = 10
search = rmsp.Search(max_comps=40)
simcache = rmsp.Simulator().simulate(
    grid.griddef, comps, 'NS_CU', search, variomodel, reals=num_reals,
    seed=515151, num_threads=8, cache=f'simcache/', cache_var='NS_CU')

Back-transform each realization to original units, applying the transform table that is stored on the nscore object.

cache = rmsp.SimCache(len(grid), 'cache/')
for ireal in range(num_reals):
    real_ns = simcache.get_real(ireal, 'NS_CU')
    real = nscore.inverse_transform(real_ns['NS_CU'])
    cache.set_real(ireal, rmsp.GeoDataFrame(real, columns=['CU']))

Visualize a couple of realizations with respect to the data.

fig, axes = rmsp.ImageGrid(2, 1, figsize=(11, 11), cbar_mode='single')
for ireal, ax in enumerate(axes):
    real = cache.get_real(ireal, 'CU', griddef=grid.griddef)
    real.sectionplot('CU', ax=ax, title=f'Realization {ireal}', **cu_pars, **grid_pars, **plot_pars)
    comps.sectionplot_draw(ax, 'CU', **cu_pars, **comps_pars)
Selected conditional realizations of copper grades.

Check CDF reproduction. The rmsp.UniStats object allows for all relevant statistics to be calculated and stored. The statistics may then be used in a variety of checks and plots without recalculation. In this case, a mean bias is shown in the realizations. This relates to stationarity issues within the simulated domain, which may be corrected through the use of trend decorrelation via TrendTransformer (a future blog topic).

sim_unistats = [rmsp.UniStats(cache.get_real(i, 'CU')['CU']) for i in range(num_reals)]
refstats = rmsp.UniStats(comps['CU'], comps['wt'])
refstats.cdfplot_checkreals(sim_unistats, title='CU')
Cumulative distribution of copper reproduction.

Interested in using RMSP for your resource modeling challenges?