Transform 2019: Integrating Striplog and GemPy

! pip install welly striplog

Authors: M. de la Varga, Evan Bianco, Brian Burnham and Dieter Werthmüller Importing GemPy

import gempy as gp

# Importing auxiliary libraries
import numpy as np
import pandas as pn
import matplotlib.pyplot as plt
import os
import welly
from welly import Location, Project
import glob
from striplog import Striplog, Legend, Decor

pn.set_option('precision', 2)

Creating striplog object

get well header coordinates

well_heads = {'alpha': {'kb_coords': (0, 0, 0)},
              'beta': {'kb_coords': (10, 10, 0)},
              'gamma': {'kb_coords': (12, 0, 0)},
              'epsilon': {'kb_coords': (20, 0, 0)}}

Reading tops file

cwd = os.getcwd()
if 'examples' not in cwd:
    data_path = os.getcwd() + '/examples'
else:
    data_path = cwd + '/..'

print(data_path+'/data/input_data/striplog_integration/*.tops')
topsfiles = glob.glob(data_path+'/data/input_data/striplog_integration/*.tops')
topsfiles

Out:

/home/travis/build/cgre-aachen/gempy/examples/integrations/../data/input_data/striplog_integration/*.tops

['/home/travis/build/cgre-aachen/gempy/examples/integrations/../data/input_data/striplog_integration/epsilon_strip.tops', '/home/travis/build/cgre-aachen/gempy/examples/integrations/../data/input_data/striplog_integration/beta_strip.tops', '/home/travis/build/cgre-aachen/gempy/examples/integrations/../data/input_data/striplog_integration/gamma_strip.tops', '/home/travis/build/cgre-aachen/gempy/examples/integrations/../data/input_data/striplog_integration/alpha_strip.tops']

Creating striplog object

my_striplogs = []

for file in topsfiles:
    with open(file) as f:
        text = f.read()
        striplog = Striplog.from_csv(text=text)
        my_striplogs.append(striplog)

striplog_dict = {'alpha': my_striplogs[1],
                 'beta': my_striplogs[2],
                 'gamma': my_striplogs[3],
                 'epsilon': my_striplogs[0]}

striplog_dict['alpha'][0]
top0.0
primary
lithoverburden
summary50.00 m of overburden
description
data
base50.0


Plot striplog

f, a = plt.subplots(ncols=4, sharey=True)

for e, log in enumerate(striplog_dict.items()):
    log[1].plot(ax=a[e], legend=None)
f.tight_layout()
plt.show()
gempy striplog

Striplog to pandas df of bottoms

rows = []
for wellname in striplog_dict.keys():
    for i, interval in enumerate(striplog_dict[wellname]):
        surface_name = interval.primary.lith
        surface_base = interval.base.middle
        x, y = well_heads[wellname]['kb_coords'][:-1]
        series = 1
        rows.append([x, y, surface_base, surface_name, series, wellname])

column_names = ['X', 'Y', 'Z', 'surface', 'series', 'wellname']
df = pn.DataFrame(rows, columns=column_names)
df
X Y Z surface series wellname
0 0 0 50.0 overburden 1 alpha
1 0 0 80.0 miguel 1 alpha
2 0 0 100.0 evan 1 alpha
3 0 0 130.0 brian 1 alpha
4 0 0 131.0 dieter 1 alpha
5 10 10 50.0 overburden 1 beta
6 10 10 75.0 miguel 1 beta
7 10 10 100.0 evan 1 beta
8 10 10 130.0 brian 1 beta
9 10 10 131.0 dieter 1 beta
10 12 0 50.0 overburden 1 gamma
11 12 0 75.0 miguel 1 gamma
12 12 0 100.0 evan 1 gamma
13 12 0 120.0 brian 1 gamma
14 12 0 121.0 dieter 1 gamma
15 20 0 50.0 overburden 1 epsilon
16 20 0 75.0 miguel 1 epsilon
17 20 0 110.0 brian 1 epsilon
18 20 0 111.0 dieter 1 epsilon


GemPy model

Create gempy model object

geo_model = gp.create_model('welly_integration')

extent = [-100, 300, -100, 200, -150, 0]
res = [60, 60, 60]

# Initializting model using the striplog df
gp.init_data(geo_model, extent, res, surface_points_df=df)

Out:

Active grids: ['regular']

welly_integration  2020-09-15 12:25
X Y Z X_r Y_r Z_r surface series id order_series smooth
0 0 0 50.0 0.44 0.47 0.25 overburden Default series 1 1 2.00e-06
1 0 0 80.0 0.44 0.47 0.44 miguel Default series 2 1 2.00e-06
2 0 0 100.0 0.44 0.47 0.56 evan Default series 3 1 2.00e-06
3 0 0 130.0 0.44 0.47 0.74 brian Default series 4 1 2.00e-06
4 0 0 131.0 0.44 0.47 0.75 dieter Default series 5 1 2.00e-06


surface series order_surfaces color id
0 overburden Default series 1 #015482 1
1 miguel Default series 2 #9f0052 2
2 evan Default series 3 #ffbe00 3
3 brian Default series 4 #728f02 4
4 dieter Default series 5 #443988 5


dec_list = []
for e, i in enumerate(striplog_dict['alpha']):
    dec_list.append(Decor({'_colour': geo_model.surfaces.df.loc[e, 'color'],
                           'width': None,
                           'component': i.primary,
                           'hatch': None}))

welly plot with gempy colors Create Decor list

dec_list = []
for e, i in enumerate(striplog_dict['alpha']):
    dec_list.append(Decor({'_colour': geo_model.surfaces.df.loc[e, 'color'],
                           'width': None,
                           'component': i.primary,
                           'hatch': None}))

# Create legend
legend = Legend(dec_list)
legend
componentcolourwidthhatch
lithoverburden
#015482NoneNone
lithmiguel
#9f0052NoneNone
lithevan
#ffbe00NoneNone
lithbrian
#728f02NoneNone
lithdieter
#443988NoneNone


Plot striplogs:

f, a = plt.subplots(ncols=4, sharey=True)

for e, log in enumerate(striplog_dict.items()):
    log[1].plot(ax=a[e], legend=legend)
f.tight_layout()
plt.show()
gempy striplog

Modifying the coordinates to make more sense

Delete points of the basement surface since we are intepolating bottoms (that surface wont exit).

Out:

True

Adding an arbitrary orientation. Remember gempy need an orientation per series

X Y Z G_x G_y G_z smooth surface
0 -500.0 1.00e-05 1.00e-05 0.0 0.0 1.0 0.01 overburden


Cell Number: mid Direction: y

Out:

<gempy.plot.visualization_2d.Plot2D object at 0x7f8564461110>
gp.set_interpolator(geo_model)

Out:

Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float64
Number of faults:  0
Compilation Done!
Kriging values:
                   values
range            5.2e+02
$C_o$            6.5e+03
drift equations      [3]

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

Out:

Lithology ids
  [5. 5. 5. ... 1. 1. 1.]
p2d = gp.plot_2d(geo_model, cell_number=[30], show_data=True, show=True)
Cell Number: 30 Direction: y gempy striplog

Out:

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

Pinch out model

As we can see the 3D model generated above does not honor the forth well lets fix it. First lets add an unconformity: between the yellow and green layer:

geo_model.add_series('Unconformity')
order_series BottomRelation isActive isFault isFinite
Default series 1 Erosion True False False
Unconformity 2 Erosion False False False


Now we set the green layer in the second series

geo_model.map_series_to_surfaces({'Uncomformity': ['brian', 'evan', 'dieter']})
surface series order_surfaces color id
0 overburden Default series 1 #015482 1
1 miguel Default series 2 #9f0052 2
2 evan Uncomformity 1 #ffbe00 3
3 brian Uncomformity 2 #728f02 4
4 dieter Uncomformity 3 #443988 5


Lastly we need to add a dummy orientation to the new series:

geo_model.add_orientations(-500, 0, -100, 'dieter', [0, 0, 1])
X Y Z G_x G_y G_z smooth surface
0 -500.0 1.00e-05 1.00e-05 0.0 0.0 1.0 0.01 overburden
1 -500.0 0.00e+00 -1.00e+02 0.0 0.0 1.0 0.01 dieter


Now we can compute:

Out:

Lithology ids
  [5. 5. 5. ... 1. 1. 1.]
p = gp.plot_2d(geo_model, cell_number=[30], show_data=True)
f, a = plt.subplots(ncols=4, sharey=True)

for e, log in enumerate(striplog_dict.items()):
    log[1].plot(ax=a[e], legend=legend)
f.tight_layout()
plt.show()
  • Cell Number: 30 Direction: y
  • gempy striplog

Getting better but not quite there yet. Since the yellow does not show up in the last well the pinch out has to happen somewhere before so lets add an artifial point to get that shape:

geo_model.add_surface_points(200, 0, -75, 'evan');
X Y Z smooth surface
0 0.0 0.0 -50.0 2.00e-06 overburden
5 100.0 100.0 -50.0 2.00e-06 overburden
10 120.0 0.0 -50.0 2.00e-06 overburden
15 200.0 0.0 -50.0 2.00e-06 overburden
1 0.0 0.0 -80.0 2.00e-06 miguel
6 100.0 100.0 -75.0 2.00e-06 miguel
11 120.0 0.0 -75.0 2.00e-06 miguel
16 200.0 0.0 -75.0 2.00e-06 miguel
2 0.0 0.0 -100.0 2.00e-06 evan
7 100.0 100.0 -100.0 2.00e-06 evan
12 120.0 0.0 -100.0 2.00e-06 evan
18 200.0 0.0 -75.0 1.00e-06 evan
3 0.0 0.0 -130.0 2.00e-06 brian
8 100.0 100.0 -130.0 2.00e-06 brian
13 120.0 0.0 -120.0 2.00e-06 brian
17 200.0 0.0 -110.0 2.00e-06 brian


gp.compute_model(geo_model)
p = gp.plot_2d(geo_model, cell_number=[30], show_data=True)
f, a = plt.subplots(ncols=4, sharey=True)

for e, log in enumerate(striplog_dict.items()):
    log[1].plot(ax=a[e], legend=legend)
f.tight_layout()
plt.show()
  • Cell Number: 30 Direction: y
  • gempy striplog

sphinx_gallery_thumbnail_number = 7

gempy striplog

Out:

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

gp.save_model(geo_model)

Total running time of the script: ( 0 minutes 21.265 seconds)

Gallery generated by Sphinx-Gallery