Note
Click here to download the full example code
Moureze¶
These two lines are necessary only if gempy is not installed
import sys, os
os.environ["THEANO_FLAGS"] = "mode=FAST_RUN,device=cpu"
# Importing gempy
import gempy as gp
# Aux imports
import numpy as np
import pandas as pn
import matplotlib.pyplot as plt
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:
Moureze_points = pn.read_csv(
'https://raw.githubusercontent.com/Loop3D/ImplicitBenchmark/master/Moureze/Moureze_Points.csv', sep=';',
names=['X', 'Y', 'Z', 'G_x', 'G_y', 'G_z', '_'], header=0, )
Sections_EW = pn.read_csv('https://raw.githubusercontent.com/Loop3D/ImplicitBenchmark/master/Moureze/Sections_EW.csv',
sep=';',
names=['X', 'Y', 'Z', 'ID', '_'], header=1).dropna()
Sections_NS = pn.read_csv('https://raw.githubusercontent.com/Loop3D/ImplicitBenchmark/master/Moureze/Sections_NS.csv',
sep=';',
names=['X', 'Y', 'Z', 'ID', '_'], header=1).dropna()
Extracting the orientatins:
mask_surfpoints = Moureze_points['G_x'] < -9999
surfpoints = Moureze_points[mask_surfpoints]
orientations = Moureze_points[~mask_surfpoints]
Giving an arbitrary value name to the surface
surfpoints['surface'] = '0'
orientations['surface'] = '0'
X | Y | Z | G_x | G_y | G_z | _ | surface | |
---|---|---|---|---|---|---|---|---|
3425 | 177.43 | 155.90 | -152.01 | -99999.0 | -99999.0 | -99999.0 | 0.05 | 0 |
3426 | 89.92 | 86.35 | -120.03 | -99999.0 | -99999.0 | -99999.0 | 0.07 | 0 |
3427 | 75.94 | 116.72 | -140.02 | -99999.0 | -99999.0 | -99999.0 | 0.78 | 0 |
3428 | 177.96 | 233.83 | -148.83 | -99999.0 | -99999.0 | -99999.0 | 0.52 | 0 |
3429 | 46.49 | 17.74 | -148.02 | -99999.0 | -99999.0 | -99999.0 | 0.11 | 0 |
X | Y | Z | G_x | G_y | G_z | _ | surface | |
---|---|---|---|---|---|---|---|---|
3409 | 47.97 | 129.89 | -132.01 | -0.45 | -0.85 | 0.27 | 0.02 | 0 |
3420 | 175.94 | 293.73 | -138.02 | 0.22 | -0.88 | 0.41 | 0.44 | 0 |
3421 | 203.97 | 367.73 | -150.01 | 0.22 | -0.97 | -0.10 | 0.71 | 0 |
3422 | 133.93 | 225.62 | -146.76 | 0.32 | -0.87 | -0.37 | 0.13 | 0 |
3430 | 290.00 | 180.00 | -103.54 | 0.13 | -0.19 | 0.97 | 0.68 | 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()
Out:
2442336
resolution_requ = [156, 206, 76]
resolution = [77, 103, 38]
resolution_low = [45, 51, 38]
geo_model = gp.create_model('Moureze')
geo_model = gp.init_data(geo_model,
extent=[-5, 305, -5, 405, -200, -50], resolution=resolution_low,
surface_points_df=surfpoints, orientations_df=orientations,
surface_name='surface',
add_basement=True)
Out:
Active grids: ['regular']
Now we can see how the data looks so far:
gp.plot_2d(geo_model, direction='y')

Out:
<gempy.plot.visualization_2d.Plot2D object at 0x7f8f6576e490>
gp.set_interpolator(geo_model,
theano_optimizer='fast_run', dtype='float64')
Out:
Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization: fast_run
Device: cpu
Precision: float64
Number of faults: 0
Compilation Done!
Kriging values:
values
range 5.4e+02
$C_o$ 6.8e+03
drift equations [3, 3]
<gempy.core.interpolator.InterpolatorModel object at 0x7f8f988b2f10>
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:
new_range = geo_model.get_additional_data().loc[('Kriging', 'range'), 'values'] * 0.2
geo_model.modify_kriging_parameters('range', new_range)
values | |
---|---|
range | 1.1e+02 |
$C_o$ | 6.8e+03 |
drift equations | [3, 3] |
gp.compute_model(geo_model, set_solutions=True, sort_surfaces=False)
Out:
Lithology ids
[2. 2. 2. ... 1. 1. 1.]
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)
gp.plot_2d(geo_model, cell_number=[16], series_n=0, show_scalar=True)

Out:
<gempy.plot.visualization_2d.Plot2D object at 0x7f8f659d7970>
gp.plot_2d(geo_model, cell_number=16, show_data=True, direction='y')

Out:
<gempy.plot.visualization_2d.Plot2D object at 0x7f8f98c835e0>
sphinx_gallery_thumbnail_number = 4

Out:
<gempy.plot.vista.GemPyToVista object at 0x7f8f6d1fc1f0>
Export data:¶
The solution is stored in a numpy array of the following shape. Axis 0 are the scalar fields of each correspondent series/faults in the following order (except basement):
order_series | BottomRelation | isActive | isFault | isFinite | |
---|---|---|---|---|---|
Default series | 1 | Erosion | True | False | False |
Basement | 2 | Erosion | False | False | False |
For the surfaces, there are two numpy arrays, one with vertices and the other with triangles. Axis 0 is each surface in the order:
surface | series | order_surfaces | color | id | |
---|---|---|---|---|---|
0 | 0 | Default series | 1 | #015482 | 1 |
1 | basement | Basement | 1 | #9f0052 | 2 |
np.save(‘Moureze_scalar’, geo_model.solutions.scalar_field_matrix) np.save(‘Moureze_ver’, geo_model.solutions.vertices) np.save(‘Moureze_edges’, geo_model.solutions.edges) gp.plot.export_to_vtk(geo_model, ‘Moureze’)
gp.save_model(geo_model)
Out:
True
Total running time of the script: ( 1 minutes 18.575 seconds)