Note
Click here to download the full example code
Coordinate system¶
Short version¶
The used coordinate system is either a
- Left-Handed System (LHS), where Easting is the \(x\)-direction, Northing the \(y\)-direction, and positive \(z\) is pointing downwards;
- Right-Handed System (RHS), where Easting is the \(x\)-direction, Northing the \(y\)-direction, and positive \(z\) is pointing upwards.
In more detail¶
The derivation on which empymod
is based ([HuTS15]) uses a right-handed
system with \(x\) to the East, \(y\) to the South, and \(z\)
downwards (ESD). In the actual original implementation of empymod
this was
changed to a left-handed system with \(x\) to the East, \(y\) to the
North, and \(z\) downwards (END). However, empymod
can equally well be
used for a coordinate system where positive \(z\) is pointing up by just
flipping \(z\), resulting in \(x\) to the East, \(y\) to the North,
and \(z\) upwards (ENU).
Left-handed system | Right-handed system | |
---|---|---|
\(x\) | Easting | Easting |
\(y\) | Northing | Northing |
\(z\) | Down | Up |
\(\theta\) | Angle E-N | Angle E-N |
\(\varphi\) | Angle down | Angle up |
There are a few other important points to keep in mind when switching between coordinate systems:
The interfaces (
depth
) have to be defined continuously increasing or decreasing, either from lowest to highest or the other way around. E.g., a simple five-layer model with the sea-surface at 0 m, a 100 m water column, and a target of 50 m 900 m below the seafloor can be defined in four ways:[0, 100, 1000, 1050]
-> LHS low to high[0, -100, -1000, -1050]
-> RHS high to low[1050, 1000, 100, 0]
-> LHS high to low[-1050, -1000, -100, 0]
-> RHS low to high
The above point affects also all model parameters (
res
,aniso
,{e;m}perm{H;V}
). E.g., for the above five-layer example this would beres = [1e12, 0.3, 1, 50, 1]
-> LHS low to highres = [1e12, 0.3, 1, 50, 1]
-> RHS high to lowres = [1, 50, 1, 0.3, 1e12]
-> LHS high to lowres = [1, 50, 1, 0.3, 1e12]
-> RHS low to high
Note that in a two-layer scenario the values are always taken as low-to-high (as it is not possible to detect the direction from only one interface).
A source or a receiver exactly on a boundary is taken as being in the lower layer. Hence, if \(z_\rm{rec} = z_0\), where \(z_0\) is the surface, then the receiver is taken as in the air in the LHS, but as in the subsurface in the RHS. Similarly, if \(z_\rm{rec} = z_\rm{seafloor}\), then the receiver is taken as in the sea in the LHS, but as in the subsurface in the RHS. This can be avoided by never placing it exactly on a boundary, but slightly (e.g., 1 mm) in the layer where you want to have it.
In
bipole
, thedip
switches sign. Correspondingly indipole
, theab
’s containing vertical directions switch the sign for each vertical component.Sign switches also occur for magnetic sources or receivers.
In this example we first create a sketch of the LHS and RHS for visualization,
followed by a few examples using dipole
and bipole
to demonstrate the
two possibilities.
import empymod
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
from matplotlib.patches import FancyArrowPatch
plt.style.use('ggplot')
RHS vs LHS¶
Comparison of the right-handed system with positive \(z\) downwards and the left-handed system with positive \(z\) upwards. Easting is always \(x\), and Northing is \(y\).
class Arrow3D(FancyArrowPatch):
"""https://stackoverflow.com/a/29188796"""
def __init__(self, xs, ys, zs):
FancyArrowPatch.__init__(
self, (0, 0), (0, 0), mutation_scale=20, lw=1.5,
arrowstyle='-|>', color='.2', zorder=100)
self._verts3d = xs, ys, zs
def draw(self, renderer):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, _ = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
FancyArrowPatch.draw(self, renderer)
def repeated(ax, pm):
"""These are all repeated for the two subplots."""
# Coordinate system
# The first three are not visible, but for the aspect ratio of the plot.
ax.plot([-2, 12], [0, 0], [0, 0], c='w')
ax.plot([0, 0], [-2, 12], [0, 0], c='w')
ax.plot([0, 0], [0, 0], [-pm*2, pm*12], c='w')
ax.add_artist(Arrow3D([-2, 14], [0, 0], [0, 0]))
ax.add_artist(Arrow3D([0, 0], [-2, 14], [0, 0]))
ax.add_artist(Arrow3D([0, 0], [0, 0], [-pm*2, pm*14]))
# Annotate it
ax.text(12, 2, 0, r'$x$')
ax.text(0, 12, 2, r'$y$')
ax.text(-2, 0, pm*12, r'$z$')
# Helper lines
ax.plot([0, 10], [0, 10], [0, 0], '--', c='.6')
ax.plot([0, 10], [0, 0], [0, -10], '--', c='.6')
ax.plot([10, 10], [0, 10], [0, 0], ':', c='.6')
ax.plot([10, 10], [10, 10], [0, -10], ':', c='.6')
ax.plot([10, 10], [0, 0], [0, -10], ':', c='.6')
ax.plot([10, 10], [0, 10], [-10, -10], ':', c='.6')
# Resulting trajectory
ax.plot([0, 10], [0, 10], [0, -10], 'C0')
# Theta
azimuth = np.linspace(np.pi/4, np.pi/2, 31)
ax.plot(np.sin(azimuth)*5, np.cos(azimuth)*5, 0, c='C5')
ax.text(3, 5, 0, r"$\theta$", color='C5', fontsize=14)
# Phi
ax.plot(np.sin(azimuth)*7, azimuth*0, -np.cos(azimuth)*7, c='C1')
ax.view_init(azim=-60, elev=20)
# Create figure
fig = plt.figure(figsize=(8, 3.5))
# Left-handed system
ax1 = fig.add_subplot(121, projection=Axes3D.name, facecolor='w')
ax1.axis('off')
plt.title('Left-handed system (LHS)\nfor positive $z$ downwards', fontsize=12)
ax1.text(7, 0, -5, r"$\varphi$", color='C1', fontsize=14)
repeated(ax1, -1)
# Right-handed system
ax2 = fig.add_subplot(122, projection='3d', facecolor='w', sharez=ax1)
ax2.axis('off')
plt.title('Right-handed system (RHS)\nfor positive $z$ upwards', fontsize=12)
ax2.text(7, 0, -5, r"$-\varphi$", color='C1', fontsize=14)
repeated(ax2, 1)
plt.tight_layout()
plt.show()
Dipole¶
A simple example using dipole
. It is a marine case with 300 meter water
depth and a 50 m thick target 700 m below the seafloor.
off = np.linspace(500, 10000, 301)
LHS¶
In the left-handed system positive \(z\) is downwards. So we have to define our model by beginning with the air layer, followed by water, background, target, and background again. This means that all our depth-values are positive, as the air-interface \(z_0\) is at 0 m.
lhs = empymod.dipole(
src=[0, 0, 100],
rec=[off, np.zeros(off.size), 200],
depth=[0, 300, 1000, 1050],
res=[1e20, 0.3, 1, 50, 1],
# depth=[1050, 1000, 300, 0], # Alternative way, LHS high to low.
# res=[1, 50, 1, 0.3, 1e20], # " " "
freqtime=1,
verb=0
)
RHS¶
In the right-handed system positive \(z\) is upwards. So we have to define our model by beginning with the background, followed by the target, background again, water, and air. This means that all our depth-values are negative.
rhs = empymod.dipole(
src=[0, 0, -100],
rec=[off, np.zeros(off.size), -200],
depth=[0, -300, -1000, -1050],
res=[1e20, 0.3, 1, 50, 1],
# depth=[-1050, -1000, -300, 0], # Alternative way, RHS low to high.
# res=[1, 50, 1, 0.3, 1e20], # " " "
freqtime=1,
verb=0
)
Compare¶
Plotting the two confirms that the results agree, no matter if we use the LHS or the RHS definition.
plt.figure(figsize=(9, 4))
ax1 = plt.subplot(121)
plt.title('Real')
plt.plot(off/1e3, lhs.real, 'C0', label='LHS +')
plt.plot(off/1e3, -lhs.real, 'C4', label='LHS -')
plt.plot(off/1e3, rhs.real, 'C1--', label='RHS +')
plt.plot(off/1e3, -rhs.real, 'C2--', label='RHS -')
plt.yscale('log')
plt.xlabel('Offset (km)')
plt.ylabel('$E_x$ (V/m)')
plt.legend()
ax2 = plt.subplot(122, sharey=ax1)
plt.title('Imaginary')
plt.plot(off/1e3, lhs.imag, 'C0', label='LHS +')
plt.plot(off/1e3, -lhs.imag, 'C4', label='LHS -')
plt.plot(off/1e3, rhs.imag, 'C1-.', label='RHS +')
plt.plot(off/1e3, -rhs.imag, 'C2:', label='RHS -')
plt.yscale('log')
plt.xlabel('Offset (km)')
plt.ylabel('$E_x$ (V/m)')
plt.legend()
ax2.yaxis.set_label_position("right")
ax2.yaxis.tick_right()
plt.tight_layout()
plt.show()
Bipole [x1, x2, y1, y2, z1, z2]¶
A time-domain example using rotated bipoles, where we define them as \([x_1, x_2, y_1, y_2, z_1, z_2]\).
times = np.linspace(0.1, 10, 301)
inp = {'freqtime': times, 'signal': 0, 'verb': 0}
lhs = empymod.bipole(
src=[-50, 50, -10, 10, 100, 110],
rec=[6000, 6100, 20, -20, 220, 200],
depth=[0, 300, 1000, 1050],
res=[1e20, 0.3, 1, 50, 1],
**inp
)
rhs = empymod.bipole(
src=[-50, 50, -10, 10, -100, -110],
rec=[6000, 6100, 20, -20, -220, -200],
depth=[0, -300, -1000, -1050],
res=[1e20, 0.3, 1, 50, 1],
**inp
)
plt.figure()
plt.plot(times, lhs, 'C0', label='LHS')
plt.plot(times, rhs, 'C1--', label='RHS')
plt.xlabel('Time (s)')
plt.ylabel('$E_x$ (V/m)')
plt.legend()
plt.show()
Bipole [x, y, z, azimuth, dip]¶
A very similar time-domain example using rotated bipoles, but this time defining them as \([x, y, z, \theta, \varphi]\). Note that \(\varphi\) has to change the sign, while \(\theta\) does not.
lhs = empymod.bipole(
src=[0, 0, 100, 10, 20],
rec=[6000, 0, 200, -5, 15],
depth=[0, 300, 1000, 1050],
res=[1e20, 0.3, 1, 50, 1],
**inp
)
rhs = empymod.bipole(
src=[0, 0, -100, 10, -20],
rec=[6000, 0, -200, -5, -15],
depth=[0, -300, -1000, -1050],
res=[1e20, 0.3, 1, 50, 1],
**inp
)
plt.figure()
plt.plot(times, lhs, 'C0', label='LHS')
plt.plot(times, rhs, 'C1--', label='RHS')
plt.xlabel('Time (s)')
plt.ylabel('$E_x$ (V/m)')
plt.legend()
plt.show()
empymod.Report()
Total running time of the script: ( 0 minutes 33.372 seconds)
Estimated memory usage: 189 MB