We here evaluate how to analyse and visualise uncertainties in a kinematic model. The basic idea is that we have a set of drillhole locations and depths and want to know how uncertain the model is at these specific locations.
The required methods are implemented in an Experiment subclass and tested here with an application to the Gippsland Basin model.
In [1]:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())
Out[1]:
In [2]:
%matplotlib inline
In [3]:
# here the usual imports. If any of the imports fails, make sure that pynoddy is installed
# properly, ideally with 'python setup.py develop' or 'python setup.py install'
import sys, os
import matplotlib.pyplot as plt
import numpy as np
# adjust some settings for matplotlib
from matplotlib import rcParams
# print rcParams
rcParams['font.size'] = 15
# determine path of repository to set paths corretly below
repo_path = os.path.realpath('../..')
import pynoddy.history
import pynoddy.experiment
rcParams.update({'font.size': 20})
In [10]:
reload(pynoddy.history)
reload(pynoddy.output)
reload(pynoddy.experiment)
# the model itself is now part of the repository, in the examples directory:
history_file = os.path.join(repo_path, "examples/GBasin_Ve1_V4_b.his")
gipps_topo_ex = pynoddy.experiment.Experiment(history = history_file)
Some basic information about the model can be obtained with:
In [11]:
print gipps_topo_ex
We can have a quick look at the model in a section view (note that Noddy is now executed in the background when required - and the output automatically generated in the required resolution):
In [12]:
gipps_topo_ex.plot_section('y')
The base plot is not very useful - but we can create a section plot with a define vertical exaggeration (keyword ve) and plot the colorbar in horizontal orientation:
In [13]:
# gipps_topo_ex.determine_model_stratigraphy()
gipps_topo_ex.plot_section('x', ve = 5, position = 'centre',
cmap = 'YlOrRd',
title = '',
colorbar = False)
gipps_topo_ex.plot_section('y', position = 100, ve = 5.,
cmap = 'YlOrRd',
title = '',
colorbar_orientation = 'horizontal')
Note: The names of the model stratigraphy (colorbar labels) are unfortunately not defined correctly in the input file - we need to fix that, then we should get useful labels, as well!
We now need to define the parameter ranges. This step can either be done through explicit definition in the notebook (see the previous notebook on the Experiment class), or a list of parameters and defined ranges plus statistics can be read in from a csv file. This enables the convenient parameter definition in a spreadsheed (for example through Excel).
In order to be read in correctly, the header should contain the labels:
In addition, it is possible to define PDF type and parameters. For now, the following settings are supported:
We can read in the parameters simply with:
In [14]:
reload(pynoddy.experiment)
# the model itself is now part of the repository, in the examples directory:
history_file = os.path.join(repo_path, "examples/GBasin_Ve1_V4_b.his")
gipps_topo_ex = pynoddy.experiment.Experiment(history = history_file)
gipps_topo_ex.load_parameter_file(os.path.join(repo_path, "examples/gipps_params.csv"))
Before generating random prerturbations, we should now store the base version so that we can always revert to it at a later stage:
In [15]:
gipps_topo_ex.freeze()
For a reproducible experiment, we can also set the random seed:
In [16]:
gipps_topo_ex.set_random_seed(12345)
And now, let's perturb the model:
In [17]:
gipps_topo_ex.random_perturbation()
Let's see what happened: we can compare the new model to the base model as we stored it before:
In [18]:
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
gipps_topo_ex.plot_section(ax = ax1, direction = 'x', model_type = "base",
colorbar = False, title = "", ve = 5.)
gipps_topo_ex.plot_section(ax = ax2, direction = 'x', colorbar = False,
title = "", ve = 5.)
In [102]:
#
# Note: keep these lines only for debugging!
#
reload(pynoddy.output)
reload(pynoddy.history)
reload(pynoddy.experiment)
# the model itself is now part of the repository, in the examples directory:
history_file = os.path.join(repo_path, "examples/GBasin_Ve1_V4_b.his")
gipps_topo_ex = pynoddy.experiment.Experiment(history = history_file)
gipps_topo_ex.load_parameter_file(os.path.join(repo_path, "examples/gipps_params.csv"))
# freeze base state
gipps_topo_ex.freeze()
# set seed
gipps_topo_ex.set_random_seed(12345)
# randomize
gipps_topo_ex.random_perturbation()
b1 = gipps_topo_ex.get_section('x', resolution = 50, model_type = 'base')
# b1.plot_section(direction = 'x', colorbar = False, title = "", ve = 5.)
b2 = gipps_topo_ex.get_section('x', resolution = 50, model_type = 'current')
In [103]:
# b1.plot_section(direction = 'x', colorbar = True, title = "", ve = 5.)
b1 -= b2
# b1.plot_section(direction = 'x', colorbar = True, title = "", ve = 5.)
print np.min(b1.block), np.max(b1.block)
In [105]:
type(b1)
Out[105]:
...and another perturbation:
In [14]:
gipps_topo_ex.random_perturbation()
In [15]:
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(311)
ax2 = fig.add_subplot(312)
ax3 = fig.add_subplot(313)
gipps_topo_ex.plot_section(ax = ax1, direction = 'x', model_type = "base",
colorbar = False, title = "", ve = 5.)
gipps_topo_ex.plot_section(ax = ax2, direction = 'x', colorbar = False,
title = "", ve = 5.)
# plot difference
In [16]:
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
gipps_topo_ex.plot_section(ax = ax1, direction = 'x', model_type = "base",
colorbar = False, title = "", ve = 5.)
gipps_topo_ex.plot_section(ax = ax2, direction = 'x', colorbar = False,
title = "", ve = 5.)
In [109]:
gipps_topo_ex.param_stats
Out[109]:
In [ ]: