/
plot_seismicevents.py
164 lines (145 loc) · 4.27 KB
/
plot_seismicevents.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
"""
Synthetic seismic
=================
This example shows how to use the :py:mod:`pylops.utils.seismicevents` module
to quickly create synthetic seismic data to be used for toy examples and tests.
"""
import matplotlib.pyplot as plt
import numpy as np
import pylops
plt.close("all")
############################################
# Let's first define the time and space axes as well as some auxiliary input
# parameters that we will use to create a Ricker wavelet
par = {
"ox": -200,
"dx": 2,
"nx": 201,
"oy": -100,
"dy": 2,
"ny": 101,
"ot": 0,
"dt": 0.004,
"nt": 501,
"f0": 20,
"nfmax": 210,
}
# Create axis
t, t2, x, y = pylops.utils.seismicevents.makeaxis(par)
# Create wavelet
wav = pylops.utils.wavelets.ricker(np.arange(41) * par["dt"], f0=par["f0"])[0]
############################################
# We want to create a 2d data with a number of crossing linear events using the
# :py:func:`pylops.utils.seismicevents.linear2d` routine.
v = 1500
t0 = [0.2, 0.7, 1.6]
theta = [40, 0, -60]
amp = [1.0, 0.6, -2.0]
mlin, mlinwav = pylops.utils.seismicevents.linear2d(x, t, v, t0, theta, amp, wav)
############################################
# We can also create a 2d data with a number of crossing parabolic events using the
# :py:func:`pylops.utils.seismicevents.parabolic2d` routine.
px = [0, 0, 0]
pxx = [1e-5, 5e-6, 1e-6]
mpar, mparwav = pylops.utils.seismicevents.parabolic2d(x, t, t0, px, pxx, amp, wav)
############################################
# And similarly we can create a 2d data with a number of crossing hyperbolic
# events using the :py:func:`pylops.utils.seismicevents.hyperbolic2d` routine.
vrms = [500, 700, 1700]
mhyp, mhypwav = pylops.utils.seismicevents.hyperbolic2d(x, t, t0, vrms, amp, wav)
############################################
# We can now visualize the different events
# sphinx_gallery_thumbnail_number = 2
fig, axs = plt.subplots(1, 3, figsize=(9, 5))
axs[0].imshow(
mlinwav.T,
aspect="auto",
interpolation="nearest",
vmin=-2,
vmax=2,
cmap="gray",
extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[0].set_title("Linear events", fontsize=12, fontweight="bold")
axs[0].set_xlabel(r"$x(m)$")
axs[0].set_ylabel(r"$t(s)$")
axs[1].imshow(
mparwav.T,
aspect="auto",
interpolation="nearest",
vmin=-2,
vmax=2,
cmap="gray",
extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[1].set_title("Parabolic events", fontsize=12, fontweight="bold")
axs[1].set_xlabel(r"$x(m)$")
axs[1].set_ylabel(r"$t(s)$")
axs[2].imshow(
mhypwav.T,
aspect="auto",
interpolation="nearest",
vmin=-2,
vmax=2,
cmap="gray",
extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[2].set_title("Hyperbolic events", fontsize=12, fontweight="bold")
axs[2].set_xlabel(r"$x(m)$")
axs[2].set_ylabel(r"$t(s)$")
plt.tight_layout()
############################################
# Let's finally repeat the same exercise in 3d
phi = [20, 0, -10]
mlin, mlinwav = pylops.utils.seismicevents.linear3d(
x, y, t, v, t0, theta, phi, amp, wav
)
fig, axs = plt.subplots(1, 2, figsize=(7, 5), sharey=True)
fig.suptitle("Linear events in 3d", fontsize=12, fontweight="bold", y=0.95)
axs[0].imshow(
mlinwav[par["ny"] // 2].T,
aspect="auto",
interpolation="nearest",
vmin=-2,
vmax=2,
cmap="gray",
extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[0].set_xlabel(r"$x(m)$")
axs[0].set_ylabel(r"$t(s)$")
axs[1].imshow(
mlinwav[:, par["nx"] // 2].T,
aspect="auto",
interpolation="nearest",
vmin=-2,
vmax=2,
cmap="gray",
extent=(y.min(), y.max(), t.max(), t.min()),
)
axs[1].set_xlabel(r"$y(m)$")
mhyp, mhypwav = pylops.utils.seismicevents.hyperbolic3d(
x, y, t, t0, vrms, vrms, amp, wav
)
fig, axs = plt.subplots(1, 2, figsize=(7, 5), sharey=True)
fig.suptitle("Hyperbolic events in 3d", fontsize=12, fontweight="bold", y=0.95)
axs[0].imshow(
mhypwav[par["ny"] // 2].T,
aspect="auto",
interpolation="nearest",
vmin=-2,
vmax=2,
cmap="gray",
extent=(x.min(), x.max(), t.max(), t.min()),
)
axs[0].set_xlabel(r"$x(m)$")
axs[0].set_ylabel(r"$t(s)$")
axs[1].imshow(
mhypwav[:, par["nx"] // 2].T,
aspect="auto",
interpolation="nearest",
vmin=-2,
vmax=2,
cmap="gray",
extent=(y.min(), y.max(), t.max(), t.min()),
)
axs[1].set_xlabel(r"$y(m)$")