<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Model set-up to simulate apparent thickness likelihood at drillhole location
Step1:  Import in Experiment object
Step2:  Adjust properties of tilt event
Step3:  Now, we can define a stochastic variable for the tilt rotation
Step4:  This method stores the current events and properties in a "base_events" attribute which is not changed in the following experiments
Step5:  Set random seed to ensure reproducibility
Step6:  Define parameter distributions for stochastic properties
Step7:  Let's set up a simple uncertainty estimation from scratch
Step8:  As a first test, we can now extract a 1-D profile
Step9:  Next step (homework for Alex)
Step10:  First step
Step11:  Extract information at the position of the sampling line
Step12:  For comparison, see the complete model section
Step13:  Homework 2 for Alex
Step14:  Load history back and check possibility to adjust thicknesses
Step15:  Test change of layer thicknesses
Step16:  Now
<ASSISTANT_TASK:>
Python Code:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())
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
repo_path = os.path.realpath('../..')
import pynoddy.history
import numpy as np
%matplotlib inline
# Combined: model generation and output vis to test:
history = "simple_model.his"
output_name = "simple_out"
#
# A general note: the 'reload' statements are only important
# for development purposes (when modules were chnaged), but not
# in required for normal execution.
#
reload(pynoddy.history)
reload(pynoddy.events)
# create pynoddy object
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 )
nm.write_history(history)
# Compute the model
reload(pynoddy)
pynoddy.compute_model(history, output_name) 
# Plot output
import pynoddy.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")
Explanation: Model set-up to simulate apparent thickness likelihood at drillhole location
End of explanation
reload(pynoddy.history)
reload(pynoddy.events)
import pynoddy.experiment
reload(pynoddy.experiment)
ex1 = pynoddy.experiment.Experiment(history)
ex1.plot_section()
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 )
tilt_options = {'name' : 'Tilt',
                'pos' : (4000, 3500, 5000),
                'rotation' : 0.,
                'plunge_direction' : 0,
                'plunge' : 0.}
nm.add_event('tilt', tilt_options)
nm.write_history(history)
# Compute the model
reload(pynoddy)
pynoddy.compute_model(history, output_name) 
# Plot output
import pynoddy.output
reload(pynoddy.output)
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('x', layer_labels = strati_options['layer_names'][::-1], 
                  colorbar = True, title="",
                  savefig = False, fig_filename = "ex01_strati.eps")
ex1 = pynoddy.experiment.Experiment(history)
ex1.plot_section()
Explanation: Import in Experiment object
End of explanation
ex1.events[2].properties['Rotation'] = 20.
ex1.plot_section('y')
Explanation: Adjust properties of tilt event:
End of explanation
ex1.freeze()
Explanation: Now, we can define a stochastic variable for the tilt rotation:
First step: "freeze" current state as base model:
End of explanation
ex1.base_events[2].properties
Explanation: This method stores the current events and properties in a "base_events" attribute which is not changed in the following experiments:
End of explanation
ex1.set_random_seed(12345)
Explanation: Set random seed to ensure reproducibility:
End of explanation
param_stats = [{'event' : 2, 
              'parameter': 'Rotation',
              'stdev': 10.0,
              'type': 'normal'}
              ]
ex1.set_parameter_statistics(param_stats)
ex1.random_draw()
ex1.events[2].properties
ex1.plot_section(colorbar=True, colorbar_orientation='horizontal')
# extract only layer 4:
l4 = ex1.get_section('y').block[:,:,:] == 4
plt.imshow(l4[:,0,:].T, origin = 'lower left', cmap = 'gray_r')
Explanation: Define parameter distributions for stochastic properties:
End of explanation
# change resolution to increase simulation speed:
resolution = 100
ex1.change_cube_size(resolution)
# initialise output variable
tmp = ex1.get_section('y')
prob_4 = np.zeros_like(tmp.block[:,:,:])
n_draws = 1000
# now: generate random models and extract blocks of layer '4'
for i in range(n_draws):
    ex1.random_draw()
    tmp = ex1.get_section('y', resolution = resolution)
    prob_4 += (tmp.block[:,:,:] == 4)
# Normalise
prob_4 = prob_4 / float(n_draws)
fig = plt.figure(figsize = (12,8))
ax = fig.add_subplot(111)
ax.imshow(prob_4.transpose()[:,0,:], 
           origin = 'lower left',
           interpolation = 'none')
plt.title("Estimated probability of unit 4")
plt.xlabel("x (E-W)")
plt.ylabel("z")
Explanation: Let's set up a simple uncertainty estimation from scratch:
End of explanation
plt.plot(prob_4[20,:,:][0], np.arange(0,50,1))
Explanation: As a first test, we can now extract a 1-D profile:
End of explanation
reload(pynoddy.experiment)
ex1 = pynoddy.experiment.Experiment(history)
Explanation: Next step (homework for Alex): extract thickness and plot as histogram:
Extract information at sampling lines
End of explanation
ex1.add_sampling_line(2500, 3500)
Explanation: First step: add sampling line at a specified (x,y) location (note: default cube size: 1 m):
End of explanation
plt.plot(ex1.get_model_lines(), np.arange(0,5000,1))
Explanation: Extract information at the position of the sampling line:
End of explanation
ex1.plot_section()
Explanation: For comparison, see the complete model section:
End of explanation
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 )
tilt_options = {'name' : 'Tilt',
                'pos' : (4000, 3500, 5000),
                'rotation' : 0.,
                'plunge_direction' : 0,
                'plunge' : 0.}
nm.add_event('tilt', tilt_options)
nm.write_history(history)
# Compute the model
reload(pynoddy)
pynoddy.compute_model(history, output_name) 
# Plot output
import pynoddy.output
reload(pynoddy.output)
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('x', layer_labels = strati_options['layer_names'][::-1], 
                  colorbar = True, title="",
                  savefig = False, fig_filename = "ex01_strati.eps")
Explanation: Homework 2 for Alex: extract value of one layer and generate statistics with dipping layers
Change layer thickness
Next step: change layer thickness in pynoddy model
End of explanation
reload(pynoddy.history)
reload(pynoddy.events)
history = 'simple_model.his'
nm = pynoddy.history.NoddyHistory(history)
history = "test.his"
nm.write_history(history)
# Compute the model
reload(pynoddy)
pynoddy.compute_model(history, output_name) 
strati.layers[-1].property_lines['Remanent Magnetization']
Explanation: Load history back and check possibility to adjust thicknesses:
End of explanation
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('x', layer_labels = strati_options['layer_names'][::-1], 
                  colorbar = True, title="",
                  savefig = False, fig_filename = "ex01_strati.eps")
reload(pynoddy.history)
reload(pynoddy.events)
history = 'simple_model.his'
nm = pynoddy.history.NoddyHistory(history)
Explanation: Test change of layer thicknesses:
End of explanation
nm.events[1].layers[5].properties['Height'] = 4200.
nm.write_history(history)
# Compute the model
reload(pynoddy)
pynoddy.compute_model(history, output_name) 
nout = pynoddy.output.NoddyOutput(output_name)
nout.plot_section('x', layer_labels = strati_options['layer_names'][::-1], 
                  colorbar = True, title="",
                  savefig = False, fig_filename = "ex01_strati.eps")
Explanation: Now: change height of one layer:
End of explanation
<END_TASK>