Chapter 6: Analyzing Topology (WIP)

import sys

import gempy as gp

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

Creating an example Model

First let’s set up a simple example model. For that we initialize the geo_data object with the correct model extent and the resolution we like. Then we load our data points from csv files and set the series and order the formations (stratigraphic pile).

# initialize geo_data object
geo_data = gp.create_data([0, 3000, 0, 20, 0, 2000], resolution=[50, 3, 67])
# import data points

X Y Z G_x G_y G_z dip azimuth polarity formation ... Y_std Z_std dip_std azimuth_std order_series isFault formation_number annotations group_id index
0 1500.000000 6.666667 990.000000 0.868243 1.000000e-07 0.496139 60.255119 90.0 1 Fault ... NaN NaN NaN NaN 1 True 1 ${\bf{x}}_{\beta \,{\bf{1}},0}$ fault NaN
1 506.333333 9.666667 1679.333333 0.258819 1.000000e-07 0.965926 15.000000 90.0 1 Layer 2 ... NaN NaN NaN NaN 2 False 2 ${\bf{x}}_{\beta \,{\bf{2}},0}$ l2_a 1.0
2 2500.000000 9.666667 911.000000 0.258819 1.000000e-07 0.965926 15.000000 90.0 1 Layer 2 ... NaN NaN NaN NaN 2 False 2 ${\bf{x}}_{\beta \,{\bf{2}},1}$ l2_a 1.0

3 rows × 22 columns

gp.set_series(geo_data, {"fault":geo_data.get_formations()[np.where(geo_data.get_formations()=="Fault")[0][0]],
                         "Rest":np.delete(geo_data.get_formations(), np.where(geo_data.get_formations()=="Fault")[0][0])},
                           order_series = ["fault", "Rest"], verbose=0, order_formations=['Fault','Layer 2', 'Layer 3', 'Layer 4', 'Layer 5'])

And quickly have a look at the data:


Then we can compile our interpolator object and compute our model:

interp_data = gp.InterpolatorData(geo_data, u_grade=[0,1])
lith_block, fault_block = gp.compute_model(interp_data)
Compiling theano function...
Compilation Done!
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float32
Number of faults:  1
gp.plot_section(geo_data, lith_block[0], 0)

Analyzing Topology

GemPy sports in-built functionality to analyze the topology of its models. All we need for this is our geo_data object, lithology block and the fault block. We input those into gp.topology_compute and get several useful outputs:

  • an adjacency graph G, representing the topological relationships of the model
  • the centroids of the all the unique topological regions in the model (x,y,z coordinates of their center)
  • a list of all the unique labels (labels_unique)
  • two look-up-tables from the lithology id’s to the node labels, and vice versa
G, centroids, labels_unique, lith_to_labels_lot, labels_to_lith_lot = gp.topology_compute(geo_data, lith_block[0], fault_block)

After computing the model topology, we can overlay the topology graph over a model section:

gp.plot_section(geo_data, lith_block[0], 0)
gp.plot_topology(geo_data, G, centroids)

So let’s say we want to check if the green layer (layer 4) is connected across the fault. For that we first need to look up which nodes belong to the layer. In this simple model we could easily do that by looking at the plot above, but we can also use the look-up-tables provided by the topology function:

dict_keys(['3', '8'])

Okay, layer 4 is represented by nodes 3 and 8. We can now put these into topology_check_adjacency function, which puts out True if the two nodes share a connection (are adjacent) and False if not:

gp.topology.check_adjacency(G, 8, 3)

We can also easily look up to which other nodes a node is adjacent:

{1: {'edge_type': 'fault'},
 2: {'edge_type': 'fault'},
 7: {'edge_type': 'stratigraphic'},
 9: {'edge_type': 'stratigraphic'}}

The adjacency dictionary of the graph shows that node 8 is connected to nodes 1, 2, 7 and 9. If we go one level deeper in the dictionary, we can access the type of connection (edge):


This way we can directly check if node 8 and 2 (or any other pair of nodes that share a connection) are connected across a fault, or just stratigraphically.