Read and Visualise Geophysical Potential-Fields

Geophysical potential fields (gravity and magnetics) can be calculated directly from the generated kinematic model. A wide range of options also exists to consider effects of geological events on the relevant rock properties. We will here use pynoddy to simply and quickly test the effect of changing geological structures on the calculated geophysical response.


In [1]:
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

In [3]:
%matplotlib inline

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


Out[2]:

Read history file from Virtual Explorer

Many Noddy models are available on the site of the Virtual Explorer in the Structural Geophysics Atlas. We will download and use one of these models here as the base model.

We start with the history file of a "Fold and Thrust Belt" setting stored on:

http://virtualexplorer.com.au/special/noddyatlas/ch3/ch3_5/his/nfold_thrust.his

The file can directly be downloaded and opened with pynoddy:


In [3]:
reload(pynoddy.history)
his = pynoddy.history.NoddyHistory(url = "http://virtualexplorer.com.au/special/noddyatlas/ch3/ch3_5/his/nfold_thrust.his")
his.determine_model_stratigraphy()


 STRATIGRAPHY
 FAULT
 FAULT
 FAULT
 FAULT
 FOLD

In [4]:
his.change_cube_size(50)

In [4]:
# Save to (local) file to compute and visualise model
history_name = "fold_thrust.his"
his.write_history(history_name)
# his = pynoddy.history.NoddyHistory(history_name)


 STRATIGRAPHY
 FAULT
 FAULT
 FAULT
 FAULT
 FOLD

In [5]:
output = "fold_thrust_out"
pynoddy.compute_model(history_name, output)

In [6]:
reload(pynoddy.output)
# load and visualise model
h_out = pynoddy.output.NoddyOutput(output)

In [8]:
# his.determine_model_stratigraphy()
h_out.plot_section('x', 
                   layer_labels = his.model_stratigraphy, 
                   colorbar_orientation = 'horizontal', 
                   colorbar=False,
                   title = '',
#                   savefig=True, fig_filename = 'fold_thrust_NS_section.eps',
                   cmap = 'YlOrRd')



In [9]:
h_out.plot_section('y', layer_labels = his.model_stratigraphy, 
                   colorbar_orientation = 'horizontal', title = '', cmap = 'YlOrRd', 
#                   savefig=True, fig_filename = 'fold_thrust_EW_section.eps',
                   ve=1.5)



In [76]:
h_out.export_to_vtk(vtk_filename = "fold_thrust")

Visualise calculated geophysical fields

The first step is to recompute the model with the generation of the geophysical responses


In [77]:
pynoddy.compute_model(history_name, output, type = 'GEOPHYSICS')

We now get two files for the caluclated fields: '.grv' for gravity, and '.mag' for the magnetic field. We can extract the information of these files for visualisation and further processing in python:


In [78]:
reload(pynoddy.output)
geophys = pynoddy.output.NoddyGeophysics(output)


240

In [78]:


In [79]:
figsize(8,8)
# imshow(geophys.grv_data, cmap = 'jet')
# define contour levels
levels = np.arange(322,344,1)
contourf(geophys.grv_data, levels, cmap = 'gray', vmin = 324, vmax = 342)
cbar = colorbar(orientation = 'horizontal')
print levels


[322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339
 340 341 342 343]

Change history and compare gravity

As a next step, we will now change aspects of the geological history (paramtereised in as parameters of the kinematic events) and calculate the effect on the gravity. Then, we will compare the changed gravity field to the original field.

Let's have a look at the properties of the defined faults in the original model:


In [80]:
for i in range(4):
    print("\nEvent %d" % (i+2))
    print "Event type:\t" + his.events[i+2].event_type
    print "Fault slip:\t%.1f" % his.events[i+2].properties['Slip']
    print "Fault dip:\t%.1f" % his.events[i+2].properties['Dip']
    print "Dip direction:\t%.1f" % his.events[i+2].properties['Dip Direction']


Event 2
Event type:	FAULT
Fault slip:	-5000.0
Fault dip:	0.0
Dip direction:	90.0

Event 3
Event type:	FAULT
Fault slip:	-3000.0
Fault dip:	0.0
Dip direction:	90.0

Event 4
Event type:	FAULT
Fault slip:	-3000.0
Fault dip:	0.0
Dip direction:	90.0

Event 5
Event type:	FAULT
Fault slip:	12000.0
Fault dip:	80.0
Dip direction:	170.0

As a simple test, we are changing the fault slip for all the faults and simply add 1000 m to all defined slips. In order to not mess up the original model, we are creating a copy of the history object first:


In [110]:
import copy
his_changed = copy.deepcopy(his)

slip_change = 2000.

his_changed.events[2].properties['Slip'] += slip_change
his_changed.events[3].properties['Slip'] += slip_change
his_changed.events[4].properties['Slip'] += slip_change
his_changed.events[5].properties['Slip'] += slip_change

We now write the adjusted history back to a new history file and then calculate the updated gravity field:


In [111]:
his_changed.write_history('fold_thrust_changed')

In [112]:
# recompute geophysical response
pynoddy.compute_model('fold_thrust_changed', 'fold_thrust_changed_out', type = 'GEOPHYSICS')

In [113]:
# load output and visualise geophysical field
geophys_changed = pynoddy.output.NoddyGeophysics('fold_thrust_changed_out')


240

In [125]:
figsize(8,8)
# imshow(geophys_changed.grv_data, cmap = 'jet')
contourf(geophys_changed.grv_data, levels, cmap = 'gray', vmin = 324, vmax = 342)
colorbar(orientation = 'horizontal')


Out[125]:
<matplotlib.colorbar.Colorbar instance at 0x113fad680>

In [115]:
figsize(8,8)
# imshow(geophys.grv_data - geophys_changed.grv_data, cmap = 'jet')
comp_levels = np.arange(-4.,4.1,0.25)
contourf(geophys.grv_data - geophys_changed.grv_data, comp_levels, cmap = 'RdBu_r')
colorbar(orientation = 'horizontal')


Out[115]:
<matplotlib.colorbar.Colorbar instance at 0x1133ddd40>

In [116]:
his.events


Out[116]:
{1: <pynoddy.events.Stratigraphy instance at 0x1063e9320>,
 2: <pynoddy.events.Fault instance at 0x10f22a200>,
 3: <pynoddy.events.Fault instance at 0x10b472950>,
 4: <pynoddy.events.Fault instance at 0x10b4727e8>,
 5: <pynoddy.events.Fault instance at 0x10b472758>,
 6: <pynoddy.events.Fold instance at 0x10b734f80>}

In [117]:
for i in range(4):
    print("Event %d" % (i+2))
    print his.events[i+2].properties['Slip']
    print his.events[i+2].properties['Dip']
    print his.events[i+2].properties['Dip Direction']


Event 2
-5000.0
0.0
90.0
Event 3
-3000.0
0.0
90.0
Event 4
-3000.0
0.0
90.0
Event 5
12000.0
80.0
170.0

In [118]:
# recompute the geology blocks for comparison:
pynoddy.compute_model('fold_thrust_changed', 'fold_thrust_changed_out')

In [119]:
geology_changed = pynoddy.output.NoddyOutput('fold_thrust_changed_out')

In [120]:
geology_changed.plot_section('x', 
                   layer_labels = his.model_stratigraphy, 
                   colorbar_orientation = 'horizontal', 
                   colorbar=False,
                   title = '',
#                   savefig=True, fig_filename = 'fold_thrust_NS_section.eps',
                   cmap = 'YlOrRd')



In [121]:
geology_changed.plot_section('y', layer_labels = his.model_stratigraphy, 
                   colorbar_orientation = 'horizontal', title = '', cmap = 'YlOrRd', 
#                   savefig=True, fig_filename = 'fold_thrust_EW_section.eps',
                   ve=1.5)



In [122]:
# Calculate block difference and export as VTK for 3-D visualisation:
import copy
diff_model = copy.deepcopy(geology_changed)
diff_model.block -= h_out.block

In [123]:
diff_model.export_to_vtk(vtk_filename = "diff_model_fold_thrust_belt")

Figure with all results


In [126]:
fig = plt.figure(figsize=(20,8))
ax1 = fig.add_subplot(131)
# original plot
levels = np.arange(322,344,1)
cf1 = ax1.contourf(geophys.grv_data, levels, cmap = 'gray', vmin = 324, vmax = 342)
# cbar1 = ax1.colorbar(cf1, orientation = 'horizontal')
fig.colorbar(cf1, orientation='horizontal')
ax1.set_title('Gravity of original model')

ax2 = fig.add_subplot(132)




cf2 = ax2.contourf(geophys_changed.grv_data, levels, cmap = 'gray', vmin = 324, vmax = 342)
ax2.set_title('Gravity of changed model')
fig.colorbar(cf2, orientation='horizontal')

ax3 = fig.add_subplot(133)


comp_levels = np.arange(-4.,4.1,0.25)
cf3 = ax3.contourf(geophys.grv_data - geophys_changed.grv_data, comp_levels, cmap = 'RdBu_r')
ax3.set_title('Gravity difference')

fig.colorbar(cf3, orientation='horizontal')

plt.savefig("grav_ori_changed_compared.eps")



In [ ]: