<jupyter_start><jupyter_text>Synthetic seismic: wedgeWe're going to make the famous wedge model, which interpreters can use to visualize the tuning effect. Then we can extend the idea to other kinds of model. Make a wedge earth model<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
length = 80  # x range
depth = 200  # z range<jupyter_output><empty_output><jupyter_text>EXERCISEMake a NumPy array of integers with these dimensions, placing a boundary at a 'depth' of 66 and another at a depth of 133.A plot of a vertical section through this array should look something like:      |      |      ---        |        |        ---          |          |<jupyter_code># YOUR CODE HERE


# We have to pass dtype=int or we get floats.
# We need ints because we're going to use for indexing later.
model = 1 + np.tri(depth, length, -depth//3, dtype=int)

plt.imshow(model)
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Now set the upper part of the model — above the wedge — to zero.<jupyter_code>model[:depth//3,:] = 0

plt.imshow(model)
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Now we can make some Vp-rho pairs (rock 0, rock 1, and rock 2).<jupyter_code>rocks = np.array([[2540, 2550],   # <-- Upper layer
                  [2400, 2450],   # <-- Wedge
                  [2650, 2800]])  # <-- Lower layer<jupyter_output><empty_output><jupyter_text>Now we can use ['fancy indexing'](http://docs.scipy.org/doc/numpy/user/basics.indexing.html) to use `model`, which is an array of 0, 1, and 2, as the indices of the rock property pairs to 'grab' from `rocks`.<jupyter_code>earth = rocks[model]<jupyter_output><empty_output><jupyter_text>Now apply `np.prod` (product) to those Vp-rho pairs to get impedance at every sample.<jupyter_code>imp = np.apply_along_axis(np.prod, arr=earth, axis=-1)<jupyter_output><empty_output><jupyter_text>Model seismic reflectionsNow 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.  EXERCISECan you write a function to compute the reflection coefficients in this model?It should implement this equation, where $Z$ is acoustic impedance and :$$ R = \frac{Z_\mathrm{lower} - Z_\mathrm{upper}}{Z_\mathrm{lower} + Z_\mathrm{upper}} $$The result should be a sparse 2D array of shape (199, 80). The upper interface of the wedge should be positive.<jupyter_code>def make_rc(imp):
    
    # YOUR CODE HERE
    
    return rc

rc = make_rc(imp)
def make_rc(imp):
    """
    Compute reflection coefficients.
    """
    upper = imp[ :-1, :]
    lower = imp[1:  , :]
    
    return (lower - upper) / (lower + upper)

rc = make_rc(imp)<jupyter_output><empty_output><jupyter_text>You should be able to plot the RC series like so:<jupyter_code>plt.figure(figsize=(8,4))
plt.imshow(rc, aspect='auto')
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>EXERCISEImplement a Ricker wavelet of frequency $f$ with amplitude $A$ at time $t$ given by:$$ \mathbf{a}(\mathbf{t}) = (1-2 \pi^2 f^2 \mathbf{t}^2) \mathrm{e}^{-\pi^2 f^2 \mathbf{t}^2} $$<jupyter_code># YOUR CODE HERE<jupyter_output><empty_output><jupyter_text>There is an implementation in `scipy.signal` but it has a 'width parameter' instead of 'frequency' so it's harder to parameterize.Instead, we'll use `bruges` to make a wavelet:<jupyter_code>from bruges.filters import ricker

f = 25  # We'll use this later.
w, t = ricker(duration=0.128, dt=0.001, f=f, return_t=True)

plt.plot(t, w)
plt.show()<jupyter_output><empty_output><jupyter_text>EXERCISEMake an RC series 200 samples long, with one positive and one negative RC. Make a corresponding time array.Pass the RC series to `np.convolve()` along with the wavelet, then plot the resulting synthetic seismogram.<jupyter_code># YOUR CODE HERE


temp = np.zeros(200)
temp[66] = 1
temp[133] = -0.5

tr = np.convolve(temp, w, mode='same')

plt.plot(tr)<jupyter_output><empty_output><jupyter_text>Synthetic wedge It's only a little trickier for us to apply 1D convolution to every trace in our 2D reflection coeeficient matrix. NumPy provides a function, `apply_along_axis()` to apply any function along any one axis of an n-dimensional array. I don't think it's much faster than looping, but I find it easier to think about.<jupyter_code>def convolve(trace, wavelet):
    return np.convolve(trace, wavelet, mode='same')

synth = np.apply_along_axis(convolve,
                            axis=0,
                            arr=rc,
                            wavelet=w)

plt.figure(figsize=(12,6))
plt.imshow(synth, cmap="Greys", aspect=0.2)
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>EXERCISEUse `ipywidgets.interact` to turn this into an interactive plot, so that we can vary the frequency of the wavelet and see the effect on the synthetic.Here's a reminder of how to use it:    from ipywidgets import interact    @interact(a=(0, 10, 1), b=(0, 100, 10))    def main(a, b):        """Do the things!"""        print(a + b)        return<jupyter_code># YOUR CODE HERE


from ipywidgets import interact

@interact(f=(4, 100, 4))
def show(f):
    w, t = ricker(duration=0.128, dt=0.001, f=f, return_t=True)
    synth = np.apply_along_axis(convolve,
                                axis=0,
                                arr=rc,
                                wavelet=w)
    plt.figure(figsize=(12,6))
    plt.imshow(synth, cmap="Greys", aspect=0.2)
    plt.colorbar()
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11391492429704397, 'fraction_numerical': 0.022530641672674837, 'mean_word_length': 3.779500430663221, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 53, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '8834', 'license': 'Apache-2.0', 'n_tokens_mistral': 1992, 'n_tokens_neox': 1822, 'n_words': 590, 'path': 'Synthetic_wedge_short.ipynb', 'repo_name': 'EvanBianco/geocomputing_demos', 'size': '5548'}