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., 03. Simulation.


The original model (Perth_Basin) hosted on 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"
fmodel = emg3d.load(data_path + fname)['model']
fgrid = fmodel.grid
Data loaded from «/home/dtr/3-GitHub/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#

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

Reproduce the model#


The following 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 = ''
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"

    url_path + path_i,
    url_path + path_o,

# 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
    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']

# Map the different series
        "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",

_ = geo_model.reorder_series(order_series)

# Drop input data from the deleted series:

# Set faults
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.

# 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


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') + topo_name, topo)
# Load the stored topography.
topo_name = 'GemPy-II-topo.npy'
topo_path = ''
topo_path += 'emg3d/external/GemPy/'+topo_name


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

# + 'GemPy-II.h5', model=fmodel)


PyVista plot#


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)

# 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
    [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)
GemPy-II model with PyVista
Wed Aug 31 21:49:34 2022 CEST
OS Linux CPU(s) 4 Machine x86_64
Architecture 64bit RAM 15.5 GiB Environment Python
File system ext4
Python 3.9.12 | packaged by conda-forge | (main, Mar 24 2022, 23:22:55) [GCC 10.3.0]
numpy 1.22.4 scipy 1.9.0 numba 0.55.2
emg3d 1.8.0 empymod 2.2.0 xarray 2022.6.0
discretize 0.8.2 h5py 3.7.0 matplotlib 3.4.3
tqdm 4.64.0 IPython 8.4.0
Intel(R) oneAPI Math Kernel Library Version 2022.0-Product Build 20211112 for Intel(R) 64 architecture applications

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

Estimated memory usage: 9 MB

Gallery generated by Sphinx-Gallery