Application 1: Two-fault model

We describe here how to generate a simple history file for computation with Noddy using the functionality of pynoddy.

pynoddy contains the functionality to generate simple models, for example to automate the model construction process, or to enable the model construction for users who are not running Windows (for which a GUI is available).

We show here the required steps to generate the model consisting of two faults that is described in the publication as application in chapter 3.1

First, some basic imports and settings for the notebook:


In [56]:
from matplotlib import rc_params

In [57]:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())


Out[57]:

In [58]:
import sys, os
import matplotlib.pyplot as plt
# adjust some settings for matplotlib
from matplotlib import rcParams
# print rcParams
rcParams['font.size'] = 15
# determine path of repository to set paths corretly below
os.chdir(r'/Users/flow/git/pynoddy/docs/notebooks/')
repo_path = os.path.realpath('../..')
import pynoddy.history

In [59]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [60]:
rcParams.update({'font.size': 20})

Defining a stratigraphy

We start with the definition of a (base) stratigraphy for the model.


In [61]:
# Combined: model generation and output vis to test:
history = "simple_model.his"
output_name = "simple_out"
reload(pynoddy.history)
reload(pynoddy.events)
nm = pynoddy.history.NoddyHistory()
# add stratigraphy
strati_options = {'num_layers' : 8,
                  'layer_names' : ['layer 1', 'layer 2', 'layer 3', 'layer 4', 'layer 5', 'layer 6', 'layer 7', 'layer 8'],
                  'layer_thickness' : [1500, 500, 500, 500, 500, 500, 500, 500]}
nm.add_event('stratigraphy', strati_options )

# write kinematic history to file:
nm.write_history(history)

We now created a kinematic history file that can be used as an input for the command-line program Noddy to calculate the resulting geological model. The functionality of the command line program is wrapped in a python function for convenience:


In [62]:
# Compute the model
reload(pynoddy)
pynoddy.compute_model(history, output_name)

Depending on the options provided in the command line program execution, several results can be obtained from the kinematic model. The most basic output provides the geological model as a 3-D discretised block model.

This geological block model can be read in with the methods of the module pynoddy.output for postprocessing and results visualisation:


In [63]:
# Plot output
reload(pynoddy.output)
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('y', layer_labels = strati_options['layer_names'][::-1], 
                  colorbar = True, title="",
                  savefig = False, fig_filename = "ex01_strati.eps")


Adding a fault event

After creating the model stratigraphy, we now add the first fault to the model.

Note: in these examples, we always start with setting up the model from the beginning for a better readability of the code. As we already created the base model with the stratigraphy, we technically wouldn't have to repeat the step again.

The required calls for a stratigraphy and a single fault are:


In [64]:
reload(pynoddy.history)
reload(pynoddy.events)
nm = pynoddy.history.NoddyHistory()
# add stratigraphy
strati_options = {'num_layers' : 8,
                  'layer_names' : ['layer 1', 'layer 2', 'layer 3', 'layer 4', 'layer 5', 'layer 6', 'layer 7', 'layer 8'],
                  'layer_thickness' : [1500, 500, 500, 500, 500, 500, 500, 500]}
nm.add_event('stratigraphy', strati_options )




# The following options define the fault geometry:
fault_options = {'name' : 'Fault_E',
                 'pos' : (6000, 0, 5000),
                 'dip_dir' : 270,
                 'dip' : 60,
                 'slip' : 1000}

nm.add_event('fault', fault_options)

Write kinematic history to file and compute the model:


In [65]:
nm.write_history(history)
# Compute the model
pynoddy.compute_model(history, output_name)

In [66]:
# Plot output
reload(pynoddy.output)
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('y', layer_labels = strati_options['layer_names'][::-1], 
                  colorbar = False, title = "",
                  savefig = False, fig_filename = "ex01_fault_E.eps")


Adding the second fault

We now add the second fault:


In [76]:
reload(pynoddy.history)
reload(pynoddy.events)
nm = pynoddy.history.NoddyHistory()
# add stratigraphy
strati_options = {'num_layers' : 8,
                  'layer_names' : ['layer 1', 'layer 2', 'layer 3', 'layer 4', 'layer 5', 'layer 6', 'layer 7', 'layer 8'],
                  'layer_thickness' : [1500, 500, 500, 500, 500, 500, 500, 500]}
nm.add_event('stratigraphy', strati_options )

# The following options define the fault geometry:
fault_options = {'name' : 'Fault_W',
                 'pos' : (4000, 3500, 5000),
                 'dip_dir' : 90,
                 'dip' : 60,
                 'slip' : 1000}

nm.add_event('fault', fault_options)

# The following options define the fault geometry:
fault_options = {'name' : 'Fault_E',
                 'pos' : (6000, 3500, 5000),
                 'dip_dir' : 270,
                 'dip' : 60,
                 'slip' : 1000}

nm.add_event('fault', fault_options)

nm.write_history(history)

In [77]:
# Compute the model
pynoddy.compute_model(history, output_name)

In [78]:
# Plot output
reload(pynoddy.output)
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('y', layer_labels = strati_options['layer_names'][::-1], 
                  colorbar = True, title="",
                  savefig = False, fig_filename = "ex01_faults_combined.eps",
                  cmap = 'YlOrRd') # note: YlOrRd colourmap should be suitable for colorblindness!


Exporting the model to VTK

The representation of models directly in Python/ IPython is (so far) limited to 2-D slices through the model model. For a model visualisation in 3-D, the model can be exported to VTK and the created VTK file can be opened in any VTK viewer (several free options are available, e.g. Paraview, Mayavi):


In [68]:
nout.export_to_vtk()