Block Model Creation#
import geolime as geo
import numpy as np
from pyproj import CRS
import pyvista as pv
pv.set_jupyter_backend('html')
geo.Project().set_crs(CRS("EPSG:20350"))
Read Surface and Drillholes files#
dh = geo.read_file("../data/domained_drillholes.geo")
domain_solid = geo.datasets.load("rocklea_dome/domain_mesh.dxf")
domain_solid.name = "channel_hg"
gs = geo.read_file("../data/surfaces.geo")
The estimation grid is created unsing created surfaces and drillholes information for elevation of origin point. The estimation grid will use 100mx100mx1m block, which is coherent with drillholes spacing. As the surfaces were interpolated on a 50mx50m grid, we just multiply the size by 2 and divide the number of cells by 2.
nu, nv = np.array(gs.shape) // 2
size_x, size_y = np.array(gs.axis()) * 2
size_z = 1
Similarly to the modelling notebook, we use the drillholes data to compute programmatically the number of cells needed to cover all the data.
z_max = np.max(dh["Z_B"])
z_min = np.min(dh["Z_E"])
z_max, z_min
(478.6, 389.2)
nw = np.ceil((z_max - z_min) / size_z)
vx_shape = [nu, nv, nw]
Create Estimation Grid#
vx = geo.Voxel(
name="Panel Grid",
shape=vx_shape,
origin=gs.origin,
axis=[
[size_x, 0, 0],
[0, size_y, 0],
[0, 0, size_z],
]
)
Similarly to Geological Modelling, Grade Modelling is more flexible if we use a BlockModel rather than a Voxel. For example, this allows to delete unnecessary block, thus reducing the memory footprint of the algorithm.
bm_panel = geo.BlockModel(name="Panel BlockModel", xyz=vx.coords())
bm_panel.element_count()
192510
Flag grid cells position compared to surfaces#
Each block can then be label depending on its geology. The geology modelling created three interface/surface separating the detritals, channel and bedrock lithology. As they are surface, each block can be flag if it is above or below a given surface. This uses the compute_block_proportion_by_gridded_surface function. The function will compute the elevation of the gridded surface at the coordinates of the block model and set it as a new property. If the parameter proportion_above is set to True, a new property will set the proportion of the block which is above this computed elevation. If the parameter is set to False, the property will indicate the proportion wich is below the computed elevation. The parameter block_discretization allows to discretise the block to compute proportion above/below the computed elevation. If the discretization is set to (1, 1, 1) the block is not discretized and only block centroids will be checked. The resulting values of the proportion property will be 0 and 1. If the discretization is set to (2, 2, 2) for example, 8 points evenly spaced inside the block will be used to compute a proportion. Note that the more discretization there is, the slower the algorithm is.
geo.compute_block_proportion_by_gridded_surface(
block_model=bm_panel,
surface=gs,
surface_z_prop="hard_top_channel",
block_elevation_name="Z_top_channel",
block_proportion_name="prop_above_channel",
proportion_above=True,
block_discretization=(3, 3, 1)
)
geo.compute_block_proportion_by_gridded_surface(
block_model=bm_panel,
surface=gs,
surface_z_prop="hard_top_bedrock",
block_elevation_name="Z_top_bedrock",
block_proportion_name="prop_above_bedrock",
proportion_above=True,
block_discretization=(3, 3, 1)
)
geo.compute_block_proportion_by_gridded_surface(
block_model=bm_panel,
surface=gs,
surface_z_prop="hard_top_detritals",
block_elevation_name="Z_top_detritals",
block_proportion_name="prop_above_detritals",
proportion_above=True,
block_discretization=(3, 3, 1)
)
Domain can then be created on the BlockModel. First numerical code are created in a property.
bm_panel["domain_code"] = 1.
bm_panel.update_property(name="domain_code", data="3", region="(Z > Z_top_bedrock) & (Z < Z_top_channel)")
bm_panel.update_property(name="domain_code", data="2", region="(Z >= Z_top_channel)")
Domain can also be flagged as region for simplicity.
bm_panel.set_region(
name="channel",
data="(Z > Z_top_bedrock) & (Z < Z_top_channel)"
)
bm_panel.set_region(
name="detritals",
data="(Z >= Z_top_channel)"
)
bm_panel.set_region(
name="bedrock",
data="(Z <= Z_top_bedrock)"
)
Now that the BlockModel is flag for the lithology, we also need to flag the points lying inside the Solid in order to only perform the estimation on those points.
Here following the same principle as compute_block_proportion_by_gridded_surface we use the compute_block_proportion_inside_solid function, which will compute for each block the proportion inside a given Solid using the same discretization method. We then create a region “channel_hg” consisting of blocks which have at least 10% inside the Solid (note that 10% was decided arbitrarily).
geo.compute_block_proportion_inside_solid(
grid_object=bm_panel,
solid=domain_solid,
block_proportion_name='in_channel_prop', # name of the new property
block_discretization=(3, 3, 3) # proportion will be calculated on 27ths
)
bm_panel.set_region(
name="channel_hg",
data="in_channel_prop >= 0.1"
)
We can now check how many block are actually inside the Solid compared to the entire BlockModel.
bm_panel.element_count(region="channel_hg")
7560
The following 3D scene will show the drillholes colored by the Iron percentage, the block model with some transparency colorcoded by the lithology and the Solid using a wireframe display also for transparency.
dh_pv = dh.to_pyvista('Fe_pct')
bm_pv = bm_panel.to_pyvista('domain_code')
domain_solid_pv = domain_solid.to_pyvista()
p = pv.Plotter()
p.add_mesh(dh_pv.tube(radius=10))
p.add_mesh(bm_pv, opacity=0.3)
p.add_mesh(domain_solid_pv, style='wireframe')
p.set_scale(zscale=10)
p.remove_scalar_bar(title='domain_code')
p.show()
Create SMU Estimation Grid#
Localising method of estimation will be used in the next sections and require BlockModel with smaller cell size, here we create a BlockModel of SMU’s being twice smaller in X and Y direction.
bm_smu = bm_panel.downscale("SMU BlockModel", discr=[2, 2, 1])
Save block models#
bm_panel.to_file("../data/block_model_panel")
bm_smu.to_file("../data/block_model_smu")