Moureze

import os

These two lines are necessary only if gempy is not installed

# Importing gempy
import gempy as gp
import gempy_viewer as gpv

# Aux imports
import numpy as np
import pandas as pd

from gempy_engine.config import AvailableBackends

Loading surface points from repository:

With pandas we can do it directly from the web and with the right args we can directly tidy the data in gempy style:

data_path = os.path.abspath('../../data/input_data/Moureze')
Moureze_points = pd.read_csv(
    filepath_or_buffer=data_path + '/Moureze_Points.csv',
    sep=';',
    names=['X', 'Y', 'Z', 'G_x', 'G_y', 'G_z', '_'],
    header=0,
)

Sections_EW = pd.read_csv(
    filepath_or_buffer=data_path + '/Sections_EW.csv',
    sep=';',
    names=['X', 'Y', 'Z', 'ID', '_'], header=1
).dropna()

Sections_NS = pd.read_csv(
    filepath_or_buffer=data_path + '/Sections_NS.csv',
    sep=';',
    names=['X', 'Y', 'Z', 'ID', '_'], header=1
).dropna()

Extracting the orientatins:

mask_surfpoints = Moureze_points['G_x'] < -9999
surface_points = Moureze_points[mask_surfpoints][::10]
orientations = Moureze_points[~mask_surfpoints][::10]

Giving an arbitrary value name to the surface

surface_points['surface'] = '0'
orientations['surface'] = '0'

surface_points.tail()

	X	Y	Z	G_x	G_y	G_z	_	surface
3370	224.000000	65.335228	-120.0000	-99999.0	-99999.0	-99999.0	0.670953	0
3386	232.000000	20.000000	-85.1014	-99999.0	-99999.0	-99999.0	0.700194	0
3396	135.939468	154.308319	-154.0210	-99999.0	-99999.0	-99999.0	0.351625	0
3411	296.000000	378.000000	-156.4050	-99999.0	-99999.0	-99999.0	0.852947	0
3424	91.960609	151.813431	-83.2593	-99999.0	-99999.0	-99999.0	0.950689	0

orientations.tail()

	X	Y	Z	G_x	G_y	G_z	_	surface
3186	139.936676	207.712677	-154.766	-0.034360	0.068784	0.997040	0.050783	0
3254	109.926086	47.649914	-135.337	0.351890	-0.766284	-0.537571	0.942560	0
3315	158.341888	211.320877	-147.312	0.053863	0.143364	0.988203	0.053091	0
3371	0.734351	90.000000	-100.000	0.756724	0.533905	0.377245	0.170960	0
3409	47.970890	129.885971	-132.010	-0.453523	-0.849487	0.269609	0.016933	0

Data initialization:

Suggested size of the axis-aligned modeling box:

Origin: -5 -5 -200

Maximum: 305 405 -50

Suggested resolution: 2m (grid size 156 x 206 x 76)

Only using one orientation because otherwhise it gets a mess

Number voxels

np.array([156, 206, 76]).prod()

2442336

resolution_requ = [156, 206, 76]
resolution = [77, 103, 38]
resolution_low = [45, 51, 38]


surface_points_table: gp.data.SurfacePointsTable = gp.data.SurfacePointsTable.from_arrays(
    x=surface_points['X'].values,
    y=surface_points['Y'].values,
    z=surface_points['Z'].values,
    names=surface_points['surface'].values.astype(str)
)

orientations_table: gp.data.OrientationsTable = gp.data.OrientationsTable.from_arrays(
    x=orientations['X'].values,
    y=orientations['Y'].values,
    z=orientations['Z'].values,
    G_x=orientations['G_x'].values,
    G_y=orientations['G_y'].values,
    G_z=orientations['G_z'].values,
    names=orientations['surface'].values.astype(str),
    name_id_map=surface_points_table.name_id_map  # ! Make sure that ids and names are shared
)

structural_frame: gp.data.StructuralFrame = gp.data.StructuralFrame.from_data_tables(
    surface_points=surface_points_table,
    orientations=orientations_table
)

geo_model: gp.data.GeoModel = gp.create_geomodel(
    project_name='Moureze',
    extent=[-5, 305, -5, 405, -200, -50],
    resolution=resolution_low,
    refinement=5,
    structural_frame=structural_frame
)

Now we can see how the data looks so far:

gpv.plot_2d(geo_model, direction='y')

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7ff62bf41cf0>

The default range is always the diagonal of the extent. Since in this model data is very close we will need to reduce the range to 5-10% of that value:

geo_model.interpolation_options.kernel_options.range *= 0.2

gp.compute_model(
    gempy_model=geo_model,
    engine_config=gp.data.GemPyEngineConfig(
        use_gpu=False,
        dtype='float32',
        backend=AvailableBackends.PYTORCH
    )
)

Setting Backend To: AvailableBackends.PYTORCH
/home/leguark/gempy/gempy/core/data/geo_model.py:164: UserWarning: You are using refinement and passing a regular grid. The resolution of the regular grid will be overwritten
  warnings.warn(
Condition number: 35643652.0.
/home/leguark/gempy_engine/gempy_engine/core/backend_tensor.py:171: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).
  return torch.tensor(array_like, dtype=dtype)

Solutions: 5 Octree Levels, 1 DualContouringMeshes

Time

300k voxels 3.5k points

• Nvidia 2080: 500 ms ± 1.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each), Memory 1 Gb

• CPU 14.2 s ± 82.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each), Memory: 1.3 Gb

2.4 M voxels, 3.5k points

• CPU 2min 33s ± 216 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) Memory: 1.3 GB

• Nvidia 2080: 1.92 s ± 6.74 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) 1 Gb

2.4 M voxels, 3.5k points 3.5 k orientations

• Nvidia 2080: 2.53 s ± 1.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

gpv.plot_2d(geo_model, cell_number='mid', series_n=0, show_scalar=True)

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7ff57ae95b10>

gpv.plot_2d(geo_model, cell_number='mid', show_data=True, direction='y')

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7ff57ae959f0>

sphinx_gallery_thumbnail_number = 4

gpv.plot_3d(geo_model)

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7ff57ae97cd0>