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'

surfpoints.tail()

	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

orientations.tail()

	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 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:

val = .1
geo_model.interpolator.theano_graph.a_T.set_value(val)
geo_model.interpolator.theano_graph.a_T_surface.set_value(val)

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 0x7ff978cff710>

gp.plot_2d(geo_model, cell_number=16, show_data=True, direction='y')

Out:

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

sphinx_gallery_thumbnail_number = 4

gp.plot_3d(geo_model)

Out:

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

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):

geo_model.series

	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:

geo_model.surfaces

	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