Hecho

These two lines are necessary only if gempy is not installed

import sys, os
import urllib

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:

dfs = []
data_amount = 'Full'

# First stratigraphic data
for letter in range(1, 10):

    try:
        dfs.append(pn.read_csv('https://raw.githubusercontent.com/Loop3D/ImplicitBenchmark/master'
                               f'/Hecho/{data_amount}/H' +
                               str(letter) + '.csv', sep=';',
                               names=['X', 'Y', 'Z', 'surface', '_'], header=0))
    except urllib.error.HTTPError as e:
        print(e, letter)

# Also faults
for f in range(1, 4):
    fault_df = pn.read_csv('https://raw.githubusercontent.com/Loop3D/ImplicitBenchmark/master/Hecho/F' +
                           str(f) + 'Line.csv',
                           sep=';', names=['X', 'Y', 'Z'], header=0)
    fault_df['surface'] = 'f' + str(f)
    dfs.append(fault_df)

# We put all the surfaces points together because is how gempy likes it:
surface_points = pn.concat(dfs, sort=True)
surface_points.reset_index(inplace=True, drop=False)
surface_points.tail()
index X Y Z _ surface
761 6 11.14 -0.17 1.53 NaN f3
762 7 11.17 0.28 1.68 NaN f3
763 8 11.16 0.40 1.82 NaN f3
764 9 11.07 0.03 1.92 NaN f3
765 10 10.89 -0.43 2.10 NaN f3


Now we do the same with the orientations:

orientations = pn.read_csv('https://raw.githubusercontent.com/Loop3D/ImplicitBenchmark/master/Hecho/Sparse/Dips.csv',
                           sep=';',
                           names=['X', 'Y', 'Z', 'G_x', 'G_z', '_'], header=0)
# Orientation needs to belong to a surface. This is mainly to categorize to which series belong and to
# use the same color
orientations['surface'] = 0

# We fill the laking direction with a dummy value:
orientations['G_y'] = 0

# Drop unecesary data point
orientations.drop([1, 3, 4], inplace=True)
orientations
X Y Z G_x G_z _ surface G_y
0 3.81 0.16 0.58 0.62 0.79 4 0 0
2 9.87 0.06 2.52 -0.76 0.64 20 0 0
5 4.83 -0.03 3.33 0.07 1.00 28 0 0


Data initialization:

Suggested size of the axis-aligned modeling box: Origin: 0 -0.5 0 Maximum: 16 0.5 4.5

Suggested resolution: 0.05m (grid size 321 x 21 x 91)

geo_model = gp.create_model('Moureze')
geo_model = gp.init_data(geo_model, extent=[0, 16, -0.5, 0.5, 0, 4.5], resolution=[321, 21, 91],
                         surface_points_df=surface_points, orientations_df=orientations,
                         surface_name='surface',
                         add_basement=True)

Out:

Active grids: ['regular']

Out:

0
X Y Z X_r Y_r Z_r G_x G_y G_z dip azimuth polarity surface series id order_series smooth
0 3.81 0.16 0.58 0.37 0.51 0.45 0.62 0.0 0.79 38.25 90.0 1 0 Default series 1 1 0.01
2 9.87 0.06 2.52 0.57 0.50 0.51 -0.76 0.0 0.64 49.86 270.0 1 0 Default series 1 1 0.01
5 4.83 -0.03 3.33 0.40 0.50 0.54 0.07 0.0 1.00 3.89 90.0 1 0 Default series 1 1 0.01


We need an orientation per series/fault. The faults does not have orientation so the easiest is to create an orientation from the surface points availablle:

f_names = ['f1', 'f2', 'f3']
for fn in f_names:
    fault_idx = geo_model.surface_points.df.index[geo_model.surface_points.df['surface'] == fn]
    gp.set_orientation_from_surface_points(geo_model, fault_idx)

Now we can see how the data looks so far:

Cell Number: mid Direction: y

Out:

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

By default all surfaces belong to one unique series.

surface series order_surfaces color id
0 0 Default series 1 #015482 1
1 0.78 Default series 2 #9f0052 2
2 1.15 Default series 3 #ffbe00 3
3 1.90 Default series 4 #728f02 4
4 2.50 Default series 5 #443988 5
5 3.10 Default series 6 #ff3f20 6
6 3.90 Default series 7 #5DA629 7
7 4.40 Default series 8 #4878d0 8
8 5.20 Default series 9 #ee854a 9
9 f1 Default series 10 #6acc64 10
10 f2 Default series 11 #d65f5f 11
11 f3 Default series 12 #956cb4 12
12 basement Basement 1 #8c613c 13


Out:

X                float64
Y                float64
Z                float64
X_r              float64
Y_r              float64
Z_r              float64
G_x              float64
G_y              float64
G_z              float64
dip              float64
azimuth          float64
polarity           int64
surface         category
series          category
id                 int64
order_series       int64
smooth           float64
dtype: object

We will need to separate with surface belong to each series:

gp.map_stack_to_surfaces(geo_model, {'Fault1': 'f1', 'Fault2': 'f2', 'Fault3': 'f3'})
surface series order_surfaces color id
0 0 Default series 1 #015482 1
1 0.78 Default series 2 #9f0052 2
2 1.15 Default series 3 #ffbe00 3
3 1.90 Default series 4 #728f02 4
4 2.50 Default series 5 #443988 5
5 3.10 Default series 6 #ff3f20 6
6 3.90 Default series 7 #5DA629 7
7 4.40 Default series 8 #4878d0 8
8 5.20 Default series 9 #ee854a 9
9 f1 Fault1 1 #6acc64 10
10 f2 Fault2 1 #d65f5f 11
11 f3 Fault3 1 #956cb4 12
12 basement Basement 1 #8c613c 13


However if we want the faults to offset the “Default series”, they will need to be more recent (higher on the pile). We can modify the order by:

geo_model.modify_order_series(4, 'Default series')
order_series BottomRelation isActive isFault isFinite
Fault3 1 Erosion True False False
Fault1 2 Erosion True False False
Fault2 3 Erosion True False False
Default series 4 Erosion True False False
Basement 5 Erosion False False False


Lastly, so far we did not specify which series/faults are actula faults:

geo_model.set_is_fault(['Fault1', 'Fault2', 'Fault3'])

Out:

Fault colors changed. If you do not like this behavior, set change_color to False.
order_series BottomRelation isActive isFault isFinite
Fault3 1 Fault True True False
Fault1 2 Fault True True False
Fault2 3 Fault True True False
Default series 4 Erosion True False False
Basement 5 Erosion False False False


Now we are good to go:

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:  3
Compilation Done!
Kriging values:
                           values
range                         17
$C_o$                        6.6
drift equations  [3, 3, 3, 3, 3]

<gempy.core.interpolator.InterpolatorModel object at 0x7ff97416ef50>

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._interpolator.theano_graph.a_T.get_value()
geo_model._interpolator.theano_graph.a_T.set_value(.2)
gp.compute_model(geo_model, sort_surfaces=True, compute_mesh=False)

Out:

Lithology ids
  [10. 10. 10. ...  8.  8.  8.]

Time

  • GTX 2080 164 ms ± 2.03 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

gp.plot_2d(geo_model, cell_number=[10], series_n=3, show_scalar=True)
Cell Number: 10 Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7ff9ae540210>
gp.plot_2d(geo_model, cell_number=[10], show_data=True)
Cell Number: 10 Direction: y

Out:

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

sphinx_gallery_thumbnail_number = 3

gp.plot_3d(geo_model, kwargs_plot_structured_grid={'opacity': 8})
Hecho

Out:

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

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
Fault3 1 Fault True True False
Fault1 2 Fault True True False
Fault2 3 Fault True True False
Default series 4 Erosion True False False
Basement 5 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
11 f3 Fault3 1 #527682 1
9 f1 Fault1 1 #527682 2
10 f2 Fault2 1 #527682 3
0 0 Default series 1 #015482 4
1 0.78 Default series 2 #9f0052 5
2 1.15 Default series 3 #ffbe00 6
3 1.90 Default series 4 #728f02 7
4 2.50 Default series 5 #443988 8
5 3.10 Default series 6 #ff3f20 9
6 3.90 Default series 7 #5DA629 10
7 4.40 Default series 8 #4878d0 11
8 5.20 Default series 9 #ee854a 12
12 basement Basement 1 #8c613c 13


def write_property_to_gocad_voxet(propertyfilename, propertyvalues):
    """
    This function writes a numpy array into the right format for a gocad
    voxet property file. This assumet there is a property already added to the .vo file,
    and is just updating the file.
    propertyfile - string giving the path to the file to write
    propertyvalues - numpy array nz,ny,nx ordering and in float format
    """
    propertyvalues = propertyvalues.astype('>f4')  # big endian
    #     array = propertyvalues.newbyteorder()
    propertyvalues.tofile(propertyfilename)
write_property_to_gocad_voxet(f'hecho_sf_gempy_{data_amount}',
                              geo_model.solutions.scalar_field_matrix[3].reshape([321, 21, 91]).ravel('F'))
4

Out:

4

Total running time of the script: ( 5 minutes 14.133 seconds)

Gallery generated by Sphinx-Gallery