Numerical examples of the book
Anton Ziolkowski and Evert Slob, 2019, Cambridge University Press; ISBN: 9781107058620.
Copyright 2018 Dieter Werthmüller
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.
import numpy as np
import empymod as epm
import matplotlib.pyplot as plt
%matplotlib inline
The book shows in the chapter black-and-white figures, and in the plate section coloured versions.
plt.rcParams['image.cmap'] = 'viridis' # Colour
# plt.rcParams['image.cmap'] = 'gray' # Black-and-white
# Offsets for period-amplitude plots
x = 5e2+np.linspace(0, 1.95e4, 196)
# Periods for offset-amplitude plots
t = np.logspace(-3, 1, 401)
# Common model parameters
model = {'src': [0, 0, 5e-50], # Source location
'rec': [x, x*0, 1e-50], # Receiver location
'freqtime': t, # Times
'signal': 0, # Impulse response
'res': [2e14, 20, 400, 20], # Resistivity
'aniso': [1, np.sqrt(2), np.sqrt(2), np.sqrt(2)], # Anisotropy
'epermH': [0, 1, 1, 1], # To avoid numerical issues due to air layer
'epermV': [0, 1, 1, 1], # " " "
'htarg': 'key_201_2012', # Use same filter as in book
# 'htarg': ['key_201_2012', -1], # Use lagged conv. Much faster.
'ft': 'sin', # Sine-Fourier
'ftarg': 'key_SinCos_201_2012', # Use same filter as in book
'verb': 1} # Verbosity
# Target at 1 km depth
ex1km = epm.dipole(depth=[0, 1000, 1040], **model)
# Target at 3 km depth
ex3km = epm.dipole(depth=[0, 3000, 3040], **model)
# Analytical, diffusive half-space solution
hs = epm.analytical(src=[0, 0, 1e-50], rec=[x, x*0, 1e-50], signal=0, res=20,
aniso=np.sqrt(2), freqtime=t, solution='dhs', verb=1)
# Generate noise
mask = 1e-2*(1.5 - np.random.random_sample(hs.shape))
mask0 = 1e-2*(1.5 - np.random.random_sample(hs.shape))
maskn = 1e-14*(1.5 - np.random.random_sample(hs.shape))
maskn0 = 1e-14*(1.5 - np.random.random_sample(hs.shape))
# Apply noise
ex1kmn = ex1km*(1 + mask) + maskn
ex3kmn = ex3km*(1 + mask) + maskn
hsn = hs*(1 + mask0) + maskn0
def plot_result(data, error=False):
"""Create figure.
Most figures in this section are very similar, and just differ in
the shown data. We can therefore define a plot-function which
we subsequently call with the different results.
"""
n = len(data)
# Define figure size depending if 2 or 4 data-sets.
if n == 2:
tit = 'Depth to resistor: '
titadd = ['1 km', '3 km']
size = (12, 6)
clbargs = {'location': 'bottom', 'fraction': .05,
'pad': 0.2, 'aspect': 30}
else:
tit = 'Half-space response'
titadd = ['', '']
size = (6, 6)
clbargs = {'location': 'bottom'}
if error:
clbticks = [-2, -1, 0, 1, 2]
vmin = -2
vmax = 2
else:
clbticks = [-16, -14, -12, -10, -8]
vmin = -16
vmax = -8
# Start figure
fig, axs = plt.subplots(figsize=size, nrows=1, ncols=n)
if n == 2:
axs = axs.ravel()
plt.subplots_adjust(hspace=0.3, wspace=0.3)
else:
axs = [axs,]
# Loop over data
for i, val in enumerate(data):
plt.sca(axs[i])
plt.title(tit + titadd[i % 2], fontweight='bold')
# pcolormesh takes x/y as start and end point of pixels,
# not midpoints. So we have to create these.
dx = (x[1]-x[0])/2
px = np.r_[x-dx, x[-1]+dx]
lt = np.log10(t)
dt = (lt[1]-lt[0])/2
pt = 10**(np.r_[lt-dt, lt[-1]+dt])
# Plot result
cs = plt.pcolormesh(px/1000, pt, np.log10(np.abs(val)),
vmin=vmin, vmax=vmax)
# Plot contours
level_factor = 2
if (vmax-vmin) < 4:
level_factor *= 2
levels = np.linspace(vmin, vmax, level_factor*(vmax-vmin)+1)
cs2 = plt.contour(x/1000, t, np.log10(np.abs(val)),
levels=levels,
linewidths=0.5, colors=[(0, 0, 0, 0.5)])
plt.xticks([5, 10, 15, 20])
plt.xlabel('Offset (km)')
plt.ylabel('$\log_{10}(t)$ (s)')
axs[i].set_yscale('log')
axs[i].invert_yaxis()
plt.yticks([1e1, 1e0, 1e-1, 1e-2, 1e-3], ('1', '0', '-1', '-2', '-3'))
# Plot colorbar
cax, kw = plt.matplotlib.colorbar.make_axes(axs, **clbargs)
cb = plt.colorbar(cs, cax=cax, ticks=clbticks, **kw)
plt.show()
def set_axis(t=False):
"""General settings for line plots."""
if np.any(t):
plt.xlabel('Time (s)')
else:
plt.xlabel('Offset (km)')
plt.xlim([0, 15])
plt.xticks([0, 5, 10, 15])
plt.ylabel('Electric field (V/m)')
plt.ylim([1e-14, plt.gca().get_ylim()[1]])
plt.legend()
plot_result([hs])
Figure 5.55 Electric field impulse response as a function of logarithmic time and offset for the VTI half-space model.
plot_result([ex1km, ex3km])
Figure 5.56 Electric field impulse response as a function of logarithmic time and offset for the land VTI model with the 40 m thick buried resistive layer at 1 km depth (left) and 3 km depth (right).
plot_result([(ex1kmn-hsn)/hsn, (ex3kmn-hsn)/hsn], True)
Figure 5.57 Electric field impulse response difference between the VTI models with and without the resistive layer as shown in Figure 5.52, but using the noise model, as a function of logarithmic time and offset for the 1 km deep (left) and 3 km deep (right) resistive layer. The colour bar indicates the normalised difference in logarithmic scale.
# Define indices
nx1, nx3 = 58, 155
nt1, nt3 = 121, 213
plt.figure(figsize=(10, 8))
plt.subplots_adjust(wspace=.3, hspace=.4)
# 1st subplot
plt.subplot(221)
plt.title('Time = '+str(np.round(t[nt1], 3))+' s; depth to resistor: 1 km',
fontweight='bold')
plt.semilogy(x/1000, ex1kmn[nt1, :], 'k-', label='With resistor')
plt.semilogy(x/1000, hsn[nt1, :], 'k--', label='Without resistor')
set_axis()
# 2nd subplot
plt.subplot(223)
plt.title('Time = '+str(np.round(t[nt3], 2))+' s; depth to resistor: 3 km',
fontweight='bold')
plt.semilogy(x/1000, ex3kmn[nt3, :], 'k-', label='With resistor')
plt.semilogy(x/1000, hsn[nt3, :], 'k--', label='Without resistor')
set_axis()
# 3rd subplot
plt.subplot(222)
plt.title('Offset = '+str(np.round(x[nx1])/1000) +
' km; depth to resistor: 1 km', fontweight='bold')
plt.loglog(t, ex1kmn[:, nx1], 'k-', label='With resistor')
plt.loglog(t, hsn[:, nx1], 'k--', label='Without resistor')
plt.xlim([4e-3, 1e1])
set_axis(t)
plt.ylim([1e-14, 2e-10])
# 4th subplot
plt.subplot(224)
plt.title('Offset = '+str(np.round(x[nx3])/1000) +
' km; depth to resistor: 3 km', fontweight='bold')
plt.loglog(t, ex3kmn[:, nx3], 'k-', label='With resistor')
plt.loglog(t, hsn[:, nx3], 'k--', label='Without resistor')
plt.xlim([1e-1, 1e1])
set_axis(t)
plt.ylim([1e-14, 1e-12])
plt.show()
Figure 5.58 Electric field impulse response for the land VTI model with (solid lines) and without (dashed lines) the resistive layer buried at 1 km (top) and 3 km (bottom) depth as a function of offset at the time of the maximum difference (left) and as a function of time at the offset of the maximum difference (right).
epm.versions('HTML')