GemPy-II: Perth Basin

This example is mainly about building a deep marine resistivity model that can be used in other examples. There is not a lot of explanation. For more details regarding the integration of GemPy and emg3d see the GemPy-I: Simple Fault Model, and make sure to consult the many useful information over at GemPy.

The model is based on the Perth Basin Model from GemPy. We take the model, assign resistivities to the lithologies, create a random topography, move it 2 km down, fill it up with sea water, and add an air layer. The result is what is referred to in other examples as model GemPy-II, a synthetic, deep-marine resistivity model.

This model is used in, e.g., 3. Simulation.

Note

The original model (Perth_Basin) hosted on https://github.com/cgre-aachen/gempy_data is released under the LGPL-3.0 License.

import os
import emg3d
import pooch
from matplotlib.colors import LogNorm

# Adjust this path to a folder of your choice.
data_path = os.path.join('..', 'download', '')

Fetch the model

Retrieve and load the pre-computed resistivity model.

fname = "GemPy-II.h5"
pooch.retrieve(
    'https://raw.github.com/emsig/data/2021-05-21/emg3d/models/'+fname,
    'ea8c23be80522d3ca8f36742c93758370df89188816f50cb4e1b2a6a3012d659',
    fname=fname,
    path=data_path,
)
fmodel = emg3d.load(data_path + fname)['model']
fgrid = fmodel.grid

Out:

Data loaded from «/home/dtr/Codes/emsig/emg3d-gallery/examples/download/GemPy-II.h5»
[emg3d v1.0.0rc3.dev5+g0cd9e09 (format 1.0) on 2021-05-21T18:40:16.721968].

QC resistivity model

fgrid.plot_3d_slicer(
    fmodel.property_x, zslice=-3000, xslice=12000,
    pcolor_opts={'norm': LogNorm(vmin=0.3, vmax=100)}
)
GemPy II

Reproduce the model

Note

The coming sections are about how to reproduce the model. For this you have to install gempy. The code example and the GemPy-II.h5-file used in the gallery were created on 2021-05-21 with gempy=2.2.9 and pandas=1.2.4.

Get and initiate the Perth Basin

import gempy as gempy
import numpy as np

# Initiate a model
geo_model = gempy.create_model('GemPy-II')
url_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/'
url_path += 'master/data/input_data/Perth_basin/'
path_i = "Paper_GU2F_sc_faults_topo_Points.csv"
path_o = "Paper_GU2F_sc_faults_topo_Foliations.csv"

pooch.retrieve(
    url_path + path_i,
    'f2964249dd941ceac35355beb78abd9c3189347fa6b845b6795240cc1a2f44d9',
    fname=path_i,
    path=data_path,
)
pooch.retrieve(
    url_path + path_o,
   'a6566d5caa8ce2fdcd4e4cb0ca643602436ca697342afacb31b7f4bd1d17c83d',
    fname=path_o,
    path=data_path,
)


# Define the grid
nx, ny, nz = 100, 100, 100
extent = [337000, 400000, 6640000, 6710000, -12000, 1000]

# Importing the data from CSV-files and setting extent and resolution
gempy.init_data(
    geo_model,
    extent=extent,
    resolution=[nx, ny, nz],
    path_i=data_path + "Paper_GU2F_sc_faults_topo_Points.csv",
    path_o=data_path + "Paper_GU2F_sc_faults_topo_Foliations.csv",
)

Initiate the stratigraphies and faults

# We just follow the example here
del_surfaces = ['Cadda', 'Woodada_Kockatea', 'Cattamarra']
geo_model.delete_surfaces(del_surfaces)

# Map the different series
gempy.map_series_to_surfaces(
    geo_model,
    {
        "fault_Abrolhos_Transfer": ["Abrolhos_Transfer"],
        "fault_Coomallo": ["Coomallo"],
        "fault_Eneabba_South": ["Eneabba_South"],
        "fault_Hypo_fault_W": ["Hypo_fault_W"],
        "fault_Hypo_fault_E": ["Hypo_fault_E"],
        "fault_Urella_North": ["Urella_North"],
        "fault_Urella_South": ["Urella_South"],
        "fault_Darling": ["Darling"],
        "Sedimentary_Series": ['Cretaceous', 'Yarragadee', 'Eneabba',
                               'Lesueur', 'Permian']
    }
)

order_series = ["fault_Abrolhos_Transfer",
                "fault_Coomallo",
                "fault_Eneabba_South",
                "fault_Hypo_fault_W",
                "fault_Hypo_fault_E",
                "fault_Urella_North",
                "fault_Darling",
                "fault_Urella_South",
                "Sedimentary_Series",
                "Basement"]

_ = geo_model.reorder_series(order_series)

# Drop input data from the deleted series:
geo_model.surface_points.df.dropna(inplace=True)
geo_model.orientations.df.dropna(inplace=True)

# Set faults
geo_model.set_is_fault(["fault_Abrolhos_Transfer",
                        "fault_Coomallo",
                        "fault_Eneabba_South",
                        "fault_Hypo_fault_W",
                        "fault_Hypo_fault_E",
                        "fault_Urella_North",
                        "fault_Darling",
                        "fault_Urella_South"])
fr = geo_model.faults.faults_relations_df.values
fr[:, :-2] = False
_ = geo_model.set_fault_relation(fr)

Compute the model with GemPy

# Set the interpolator.
gempy.set_interpolator(
    geo_model,
    compile_theano=True,
    theano_optimizer='fast_run',
    gradient=False,
    dtype='float32',
    verbose=[]
)

# Compute it.
sol = gempy.compute_model(geo_model, compute_mesh=True)

# Get the solution at the internal grid points.
sol = gempy.compute_model(geo_model)

Assign resistivities to the id’s

We define here a discretize mesh identical to the mesh used by GemPy, and subsequently assign resistivities to the different lithologies.

Please note that these resistivities are made up, and do not necessarily relate to the actual lithologies.

# We create a mesh 20 km x 20 km x 5 km, starting at the origin.
# As long as we have the same number of cells we can trick the grid
# original into any grid we want.
hx = np.ones(nx)*20000/nx
hy = np.ones(ny)*20000/ny
hz = np.ones(nz)*5000/nz
grid = emg3d.TensorMesh([hx, hy, hz], origin=(0, 0, -5000))

# Make up some resistivities that might be interesting to model.
ids = np.round(sol.lith_block)
res = np.ones(grid.n_cells)
res[ids == 9] = 2.0    # Cretaceous
res[ids == 10] = 1.0   # Yarragadee
res[ids == 11] = 4.0   # Eneabba
res[ids == 12] = 50.0  # Lesueur
res[ids == 13] = 7.0   # Permian
res[ids == 14] = 10.0  # Basement

Topography

Calls to geo_model.set_topography(source='random') create a random topography every time. In order to have it reproducible we saved it once and load it now.

Originally it was created and stored like this:

out = geo_model.set_topography(source='random')
np.save(data_path + topo_name, topo)
# Load the stored topography.
topo_name = 'GemPy-II-topo.npy'
topo_path = 'https://raw.github.com/emsig/data/2021-05-21/'
topo_path += 'emg3d/external/GemPy/'+topo_name

pooch.retrieve(
    topo_path,
    '10bb3d672ba26f6d8cb85eb33086daebb1c19bcbf9547c0b17d93f1c0dcf4e20',
    fname=topo_name,
    path=data_path,
)

out = geo_model.set_topography(
        source='saved', filepath=data_path+topo_name, allow_pickle=True)
topo = out.topography.values_2d

# Apply the topography to our resistivity cube.
res = res.reshape(grid.shape_cells, order='C')

# Get the scaling factor betw. original extent and our made-up extent.
fact = 5000/np.diff(extent[4:])

# Loop over all x-y-values and convert cells above topography to water.
for ix in range(nx):
    for iy in range(ny):
        res[ix, iy, grid.cell_centers_z > topo[ix, iy, 2]*fact] = 0.3

Extend the model by sea water and air

# Create an emg3d-model.
model = emg3d.Model(grid, property_x=res.ravel('F'))

# Add 2 km water and 500 m air.
fhz = np.r_[np.ones(nz)*5000/nz, 2000, 500]
z0 = -7000

# Make the full mesh
fgrid = emg3d.TensorMesh([hx, hy, fhz], origin=(0, 0, z0))

# Extend the model.
fmodel = emg3d.Model(fgrid, np.ones(fgrid.shape_cells))
fmodel.property_x[:, :, :-2] = model.property_x
fmodel.property_x[:, :, -2] = 0.3
fmodel.property_x[:, :, -1] = 1e8

# emg3d.save(data_path + 'GemPy-II.h5', model=fmodel)

fgrid

PyVista plot

Note

The final cell is about how to plot the model in 3D using PyVista, for which you have to install pyvista.

The code example was created on 2021-05-21 with pyvista=0.30.1.

import pyvista
import numpy as np

dataset = fgrid.toVTK({'res': np.log10(fmodel.property_x.ravel('F'))})

# Create the rendering scene and add a grid axes
p = pyvista.Plotter(notebook=True)
p.show_grid(location='outer')

# Add spatially referenced data to the scene
dparams = {'rng': np.log10([0.3, 500]), 'show_edges': False}
xyz = (17500, 17500, -1500)
p.add_mesh(dataset.slice('x', xyz), name='x-slice', **dparams)
p.add_mesh(dataset.slice('y', xyz), name='y-slice', **dparams)

# Add a layer as 3D
p.add_mesh(dataset.threshold(
    [np.log10(49.9), np.log10(50.1)]), name='vol', **dparams)

# Show the scene!
p.camera_position = [
      (-10000, -41000, 8500), (10000, 10000, -3250), (0, 0, 1)
]
p.show()
GemPy-II model with PyVista
Sat Jul 17 18:22:01 2021 CEST
OS Linux CPU(s) 4 Machine x86_64
Architecture 64bit RAM 15.5 GiB Environment Python
Python 3.9.4 | packaged by conda-forge | (default, May 10 2021, 22:13:33) [GCC 9.3.0]
numpy 1.21.0 scipy 1.7.0 numba 0.53.1
emg3d 1.1.0 empymod 2.1.2 xarray 0.18.2
discretize 0.7.0 h5py 3.3.0 matplotlib 3.4.2
tqdm 4.61.2
Intel(R) oneAPI Math Kernel Library Version 2021.2-Product Build 20210312 for Intel(R) 64 architecture applications


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

Estimated memory usage: 25 MB

Gallery generated by Sphinx-Gallery