<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 EXERCISE — Synthetic seismic
Step1:  Now set the upper part of the model — above the wedge — to zero (blue).
Step2:  Now we can make some Vp-rho pairs (rock 0 and rock 1).
Step3:  I was using np.take here, but 'fancy indexing' is shorter and more intuitive. We are just going to index rocks using the integers in model. That is, if model has a 1, we take the second element, [2400, 2450], from rocks. We'll end up with an array containing the rocks corresponding to each element of earth.
Step4:  Now apply np.product to those Vp-rho pairs to get impedance at every sample.
Step5:  Model seismic reflections
Now we have an earth model — giving us acoustic impedance everywhere in this 2D grid — we define a function to compute reflection coefficients for every trace.
Step6:  We'll use scipy.signal.ricker to make a wavelet
Step7:  Now apply 1D convolution to every trace
Step8:  Measure amplitude at the peak
Step9:  Find the tuning thickness
We can find the thickness at which the amplitude peak starts... this is called the tuning thickness.
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt
import numpy as np
% matplotlib inline
import matplotlib as mpl
img = mpl.image.imread('../data/Hubbard_etal_2014.png')
meanimg = np.mean(img, axis=-1)
model = np.digitize(meanimg, bins=[0.5, 0.6, 0.7, 0.8])
np.unique(model)
length = 80  # x range
depth = 200  # z range
model = 1 + np.tri(depth, length, -depth//3, dtype=int)
plt.imshow(model)
plt.show()
model[100]
Explanation: EXERCISE — Synthetic seismic: wedge
Q. Try changing the rock properties.
Q. Try plotting unresolved and resolved traces.
Q. try using a different wavelet, e.g. spike, square, sinc, Ormsby, Butterworth, etc.
Q. Try making a different geometry.
Make an earth model
End of explanation
model[:depth//3,:] = 0
plt.imshow(model)
plt.colorbar()
plt.show()
Explanation: Now set the upper part of the model — above the wedge — to zero (blue).
End of explanation
rocks = np.array([[2600, 2700],
                  [2700, 2750],
                  [2400, 2450],
                  [2550, 2650],
                  [2800, 3000]])
Explanation: Now we can make some Vp-rho pairs (rock 0 and rock 1).
End of explanation
earth = rocks[model]
Explanation: I was using np.take here, but 'fancy indexing' is shorter and more intuitive. We are just going to index rocks using the integers in model. That is, if model has a 1, we take the second element, [2400, 2450], from rocks. We'll end up with an array containing the rocks corresponding to each element of earth.
End of explanation
imp = np.apply_along_axis(np.product, -1, earth)
plt.imshow(imp)
plt.colorbar()
Explanation: Now apply np.product to those Vp-rho pairs to get impedance at every sample.
End of explanation
def make_rc(imp):
    upper = imp[ :-1, :]
    lower = imp[1:  , :]
    
    return (lower - upper) / (lower + upper)
rc = make_rc(imp)
plt.imshow(rc, aspect=1)
plt.show()
Explanation: Model seismic reflections
Now we have an earth model — giving us acoustic impedance everywhere in this 2D grid — we define a function to compute reflection coefficients for every trace.
End of explanation
from scipy.signal import ricker
f = 25
wavelet = ricker(100, a=5)
plt.plot(wavelet)
plt.show()
import bruges
w = bruges.filters.ricker(duration=0.1, dt=0.001, f=60)
plt.plot(w)
Explanation: We'll use scipy.signal.ricker to make a wavelet:
End of explanation
def convolve(trace):
    return np.convolve(trace, w, mode='same')
synth = np.apply_along_axis(convolve,
                            axis=0,
                            arr=rc)
plt.figure(figsize=(12,6))
plt.imshow(synth, cmap="RdBu", aspect=1)
plt.colorbar()
plt.show()
Explanation: Now apply 1D convolution to every trace:
End of explanation
apparent = np.amin(synth, axis=0)
t_apparent = np.argmin(synth, axis=0)
actual = synth[depth//3, :]
t_actual = np.ones_like(t_apparent) * depth//3
plt.plot(apparent, label="apparent")
plt.plot(actual, label="actual")
plt.text(6, 0, "actual", color='g')
plt.text(1, -0.18, "apparent", color='b')
# plt.legend()
plt.show()
# t_apparent = t_apparent.astype(np.float)
# t_apparent[0] = np.nan
plt.figure(figsize=(12,6))
# Plot seismic
plt.imshow(synth, cmap="Greys", aspect=0.2)
plt.colorbar()
# Plot horizons
plt.plot(t_apparent, 'r')
plt.plot(t_actual, 'c')
plt.xlim([0, 79])
plt.ylim([199, 0])
plt.title("Wedge model, {} Hz wavelet".format(f))
plt.show()
Explanation: Measure amplitude at the peak
End of explanation
tuning = np.argmin(actual)
vp = 2400
wavelength = 2 * 1000 * (vp / f) / vp  # in milliseconds TWT
print("Tuning thickness is {:.1f} ms, or 1/{:.1f} wavelength".format(tuning, wavelength/tuning))
Explanation: Find the tuning thickness
We can find the thickness at which the amplitude peak starts... this is called the tuning thickness.
End of explanation
<END_TASK>