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:

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')
Cell Number: mid Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7ff9ba4f1f50>
gp.set_interpolator(geo_model,
                    theano_optimizer='fast_run', dtype='float32')

Out:

Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_run
Device:  cpu
Precision:  float32
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 0x7ff9bbbecb50>

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:

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)
Cell Number: 16 Direction: y

Out:

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

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7ff977642390>

sphinx_gallery_thumbnail_number = 4

Moureze

Out:

<gempy.plot.vista.GemPyToVista object at 0x7ff97845db90>

image0

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: ( 4 minutes 37.580 seconds)

Gallery generated by Sphinx-Gallery