Geological Modelling

Geological Modelling#

This notebook documents several workflows for modelling and domaining of data.

We will model the bulk geometry of our paleochannel using.

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

pv.set_jupyter_backend('html')
geo.Project().set_crs(CRS("EPSG:20350"))

We here create a function to wrap code that was detailed Notebook 2, which intends to project all data on a given axis (X, Y, or Z). The flatten parameter allows to project the data according to real depth or relative depth from Collar.

def project_section(
    dh: geo.Drillholes,
    coord: geo.Coord,
    var: str,
    region: str = None,
    flatten = False
):
    Z = "Z_M"
    yaxis = True

    if flatten:
        Z = "FROM"
        yaxis = "reversed"
    
    f = geo.scatter(geo_object=dh, property_x=f"{coord}_M", property_y=Z, region=region, yaxis_autorange=yaxis)
    
    # color options
    f["data"][0]["marker"]["color"] = dh.data(var, region)
    f["data"][0]["marker"]["colorscale"] = "rainbow"
    f["data"][0]["marker"]["line"]["width"] = 0
    
    return f

GeoLime drillholes file has been exported in the previous notebook. The read_file function reads any .geo. file and will load the object directly.

We reload drillhole data from previous notebook.

dh = geo.read_file("../data/dh_hyper.geo")

Create Channel Geometry#

Attention

Lithology were modelled using information provided in article from Haest, 2012.

The objective of the following section is to use geometrical shapes to model the paelo-channel based. Though this method should not be used for a resource model used in reporting, it allows us to add lithology information to the data in order to display the capabilities of the surface modelling algorithm.

project_section(dh=dh, coord="X", var="Fe_pct", flatten=True)

In order to better represent the channel geometry we use two ellipsoids allowing us to characterise the bending of the channel. Using information provided in the article the ellispoid parameters were chosen to match as best as possible the real channel geometry as we do not possess the dxf file.

First the two ellipsoid are defined by their center, here X and Y locations have been determined visually from the article.

cx = 547895.7
cy = 7475196.5
cz = np.max(dh["Z_B"])

cx_1 = 547200
cy_1 = 7473400
cz_1 = np.max(dh["Z_B"])

The the ellipsoid semi axes are also determined given the visual shapes from the image in the article.

rx = 1000
ry = 2400
rz = 75

rx_1 = 1200
ry_1 = 1200
rz_1 = 75

The first ellipsoid is an ellipsoid oriented North-South, and the second has an angle of N060, allowing it to capture channel bending.

a_yx = np.radians(120)
sina_yx =np.sin(a_yx)
cosa_yx = np.cos(a_yx)

The following code create two properties of distance based on the ellipsoids and their centroids. Each drillholes point inside the ellipsoids will have a distance less than 1. Distances computation is based on the ellipsoid formula :

\( \frac{{\left( {x - x_0 } \right)^2 }}{{a^2 }} + \frac{{\left( {y - y_0 } \right)^2 }}{{b^2 }} = 1 \)

Distance are multiply by a negative factor, in order to only keep values higher than -1, being the values inside the ellipsoids.

dh.set_property(
    name="ellipsoidal_distance_1",
    data=f"-1. * (((X_M - {cx_1})*{cosa_yx} + (Y_M - {cy_1})*{sina_yx})**2 / ({rx_1} * {rx_1}) + ((X_M - {cx_1})*{sina_yx} + (Y_M - {cy_1})*{cosa_yx} )**2 / ({ry_1} * {ry_1}) + (Z_M - {cz_1} )**2/ ({rz_1} * {rz_1}))"
)

dh.set_property(
    name="ellipsoidal_distance_2",
    data=f"-1. * ((X_M - {cx}) * (X_M - {cx}) / ({rx} * {rx}) + (Y_M - {cy}) * (Y_M - {cy}) / ({ry} * {ry}) + (Z_M - {cz}) * (Z_M - {cz}) / ({rz} * {rz}))"
)

The channel property is then created using the points located inside the two ellipsoids.

dh.set_property(
    name="channel",
    data="(((ellipsoidal_distance_1 > -1) * (-ellipsoidal_distance_1)) + ((ellipsoidal_distance_2 > -1) * (-ellipsoidal_distance_2)))/1.7 "
)

Channel shape can be visualised in 2D. Here let’s see the mean grade per drillholes where at least one sample is located within the channel.

dh.plot_2d(property="Fe_pct", agg_method="mean", region="channel > 0", region_agg="any", interactive_map=True, style_kwds={"fillColor": "none"})

Make this Notebook Trusted to load map: File -> Trust Notebook

Using the project section function the shape can also be visualised vertically along specific axis.

project_section(dh=dh, coord="X", var="channel", flatten=True)

Create Domain based on distance to channel and slope#

In order to add more complexity to the dataset, points outside the channel are also split into two categories:

detritals, being the overburden made of loose materials
bedrock

These two denomination are taken from the article.

In orde to facilitate computation and visualisation this lithology are stored as noun and as category with numeric values.

Again, due to the lack of data, the split between detritals and bedrock was made to match visually graphics from the article.

dh["domain_code"] = 1.
dh.update_property(name="domain_code", data="3.", region="channel > 0")
dh.update_property(name="domain_code", data="2.", region="(channel > 0) & (FROM < 0.25 + 7 *(7479595-Y_M)/(7479595-7472800))")

project_section(dh=dh, coord="X", var="domain_code")

project_section(dh=dh, coord="Y", var="domain_code")

Domain code are translated to litteral character

dh["domain"] = "'channel'"
dh.update_property(name="domain", data="'bedrock'", region="domain_code == 1")
dh.update_property(name="domain", data="'detritals'", region="domain_code == 2")

Analyse chemicals trend using contact analysis#

geo.contact_analysis(dh, 'domain', 'detritals', 'channel', 'Fe')

geo.contact_analysis(dh, 'domain', 'channel', 'bedrock', 'Fe')

geo.contact_analysis(dh, 'domain', 'detritals', 'channel', 'SiO2')

geo.contact_analysis(dh, 'domain', 'channel', 'bedrock', 'SiO2')

Surfaces Creation#

Now that we have lithology information on the drillholes, we intend to interpolate the surface markers onto a regulard grid. This will be achieved using the ClosestPointValue methodology ensuring to avoid pinch-out and clipp-in of shorter drillholes.

Creation of Gridded Surface for interpolation#

The grid used for the interpolation is here created using the bouding box of the drillholes data, ensuring the validity of the grid extent.

To do so we first initialise a 2D Voxel object. Voxel can be seen as shoebox grid where all cells are defined. This allows faster computation for algorithms but can lead to bigger memory footprint of object.

The Voxel is then converted into a GriddedSurface object, which can ben seen as 2D object where the elevation is a property. For example, LiDAR data can be loaded as GriddedSurface as for each coordinates (X,Y) there is only one Z value.

The 2D Voxel is created with a cell-size of 50m. This size was chosen regarding the size of the future BlockModel of Estimation: there is no need to have a better resolution than the object to filter.

size_x = 50
size_y = 50

The number of cells (shape) and the origin of the Voxel are computed automatically using the drillholes information and the provided cell size. For sanity reason here offset the origin in order to avoid having drillholes data on the edges of the Voxel, and we add 4 cells in X and Y direction.

First the bounds of Drillholes are computed. Bounds can be accessed for any GeoLime object and represent the 2 points of the minimum and maximum values of each coordinates.

dh_bbox = dh.bounds()
dh_bbox

array([[5.4549490e+05, 7.4727934e+06, 3.8920000e+02],
       [5.4849810e+05, 7.4795953e+06, 4.7860000e+02]])

The shape of the Voxel is computed based on the distances between the bounds in each direction the manually given size. Dividing the distance by the size will give the number of cell in each direction.

vx_shape = (dh_bbox[1] - dh_bbox[0]) / np.array([size_x, size_y, 1])
vx_shape

array([ 60.064, 136.038,  89.4  ])

However the shape might not be an actual integer due to the different paramters. Using numpy and function such as round/ceil/floor will transforms the floating array into integer array. Using the ceil function will take the closest upper integer, hence making sure the Voxel covers the entirety of the Drillholes

vx_shape = np.ceil(vx_shape)
vx_shape

array([ 61., 137.,  90.])

Here for gelogical modelling we only need a 2D Voxel, so we take only the two first values of the array.

vx_shape = vx_shape[:2]

For sanity reason we might also add some cells in each direction in order ensure Drillholes points not lying too close to the edge of the Voxel.

vx_shape = vx_shape + 2
vx_shape

array([ 63., 139.])

We can do the same sanity operation for the origin of the Voxel and move to about two cells from the center of the Drillholes. We also use the floor method to use the closest smaller integer.

vx_origin = np.floor(dh_bbox[0] - np.array([100, 100, 0]))
vx_origin

array([5.453940e+05, 7.472693e+06, 3.890000e+02])

We can then create a 2D Voxel using the Voxel object.

vx = geo.Voxel(
    name="Grid2D",
    shape=vx_shape,
    origin=vx_origin,
    axis=[
        [size_x, 0, 0],
        [0, size_y, 0]
    ]
)

For Geological Modelling, the Voxel has to be converted to a GriddedSurface. GriddedSurface are more suited for topography/interface modelling.

gs = geo.GriddedSurface(name="geol_surfaces", xyz=vx.coords())
print(f"Gridded Surface Shape: {gs.shape}")

Gridded Surface Shape: [ 63 139]

Creation of lithological sequence#

In order to use the ClosestPointValue methodology we need to provide the stratigraphic sequence or the lithology relationship of the different values from top to depth.

lithos = ["detritals", "channel", "bedrock"]

Estimate Surface#

The ClosestPointValue methodology first compute the constraints given the drillholes lithology values. Constraints ensure conformity of lithologies.

Drillholes missing lithology values are interpolated using RBF and must respect existing constraints. For better accuracy lithology are interpolated using thickness instead of elevation allowing to avoid issue related to trend or periodic variations.

estimator = geo.ClosestPointValue(
    method=geo.SurfaceEstimationMethod.THICKNESS,
    interpolator=geo.DrillholesInterpolationMethod.RBF
)
estimator.build_constraints(data=dh, domain_col="domain", domain_sequence=lithos)
estimator.fill_constraints()

Once the missing drillholes values have been interpolated, grid values are also interpolated using again RBF method. RBF allows to have smooth surface. Once interpolated, values ares again verified in order to respect conformity relationship.

estimator.surface_creation(grid_surf=gs, interpolator=geo.ExternalSurfaceInterpolationMethod.RBF)

Visual check of created surfaces#

gs.plot_2d(property="hard_thick_channel", agg_method="max", interactive_map=True, interactive_max_sample=9000)