<jupyter_start><jupyter_text>X lines of Python Wedge modelThis is part of [an Agile blog series](http://ageo.co/xlines00) called **x lines of Python**.We start with the usual preliminaries.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>Make an earth modelWe'll start off with an earth model --- an array of 'cells', each of which has some rock properties.Line 1 sets up some basic variables, then in line 2 I've used a little matrix-forming trick, `np.tri(m, n, k)`, which creates an *m* &times; *n* matrix with ones below the *k*th diagonal, and zeros above it. The `dtype` specification just makes sure we end up with integers, which we need later for the indexing trick.Then line 3 just sets every row above `depth//3` (the `//` is integer division, because NumPy prefers integers for indexing arrays), to 0.<jupyter_code>length, depth = 40, 100
model = 1 + np.tri(depth, length, -depth//3, dtype=int)
model[:depth//3,:] = 0<jupyter_output><empty_output><jupyter_text>We'll have a quick look with some very basic plotting commands.<jupyter_code>plt.imshow(model, cmap='viridis', aspect=0.2)
plt.show()
model[60]<jupyter_output><empty_output><jupyter_text>Now we can make some Vp-rho pairs (rock 0, rock 1, rock 2) and select from those with `np.take`. This works like `vlookup` in Excel --- it says "read this array, `model` in this case, in which the values *i* are like 0, 1, ... n, and give me the *i*th element from this other array, `rocks` in this case.<jupyter_code>rocks = np.array([[2700, 2750],  # Vp, rho
                  [2400, 2450],
                  [2800, 3000]])<jupyter_output><empty_output><jupyter_text>**Edit:** I was using `np.take` here, but ['fancy indexing'](http://docs.scipy.org/doc/numpy/user/basics.indexing.html) 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`.<jupyter_code>earth = rocks[model]<jupyter_output><empty_output><jupyter_text>Now apply `np.product` to those Vp-rho pairs to get impedance at every sample.This might look a bit magical, but we're just telling Python to apply the function `product()` to every set of numbers it encounters on the last axis (index `-1`). The array `earth` has shape (100, 40, 2), so you can think of it as a 100 row x 40 column 'section' in which each 'sample' is occupied by a Vp-rho pair. That pair is in the last axis. So product, which just takes a bunch of numbers and multiplies them, will return the impedance (the product of Vp and rho) at each sample location. We'll end up with a new 100 x 40 'section' with impedance at every sample.<jupyter_code>imp = np.apply_along_axis(np.product, -1, earth)<jupyter_output><empty_output><jupyter_text>We could have saved a step by taking from `np.product(rocks, axis=1)` but I like the elegance of having an earth model with a set of rock properties at each sample location. That's how I think about the earth --- and it's similar to the concept of a geocellular model.  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.I love this indexing trick though I admit it looks weird the first time you see it. It's easier to appreciate for a 1D array. Let's look at the differences:    >>> a = np.array([1,1,1,2,2,2,3,3,3])    >>> a[1:] - a[:-1]    array([0, 0, 1, 0, 0, 1, 0, 0])This is equivalent to:    >>> np.diff(a, axis=0)    But I prefer to spell it out so it's analogous to the sum on the denominator.<jupyter_code>rc =  (imp[1:,:] - imp[:-1,:]) / (imp[1:,:] + imp[:-1,:])

plt.imshow(rc, cmap='Greys', aspect=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>We'll use a wavelet function from [`bruges`](https://github.com/agile-geoscience/bruges). This is not cheating! Well, I don't think it is... we could use `scipy.signal.ricker` but I can't figure out how to convert frequency into the 'width' parameter that function wants. Using the Ricker from `bruges` keeps things a bit simpler.<jupyter_code>import bruges

w = bruges.filters.ricker(duration=0.100, dt=0.001, f=40)<jupyter_output><empty_output><jupyter_text>Let's make sure it looks OK:<jupyter_code>plt.plot(w)
plt.show()<jupyter_output><empty_output><jupyter_text>Now one more application of `apply_along_axis`. We could use a loop to step over the traces, but the rule of thumb in Python is "if you are using a loop, you're doing it wrong.". So, we'll use `apply_along_axis`.It looks a bit more complicated this time, because we can't just pass a function like we did with `product` before. We want to pass in some more things, not just the trace that `apply_along_axis` is going to send it. So we use Python's 'unnamed function creator', `lambda` (in keeping with all things called `lambda`, it's a bad name that no-one can quite explain).<jupyter_code>synth = np.apply_along_axis(lambda t: np.convolve(t, w, mode='same'),
                            axis=0,
                            arr=rc)

plt.imshow(synth, cmap="Greys", aspect=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>That's it! And it only needed 9 lines of Python! Not incldung boring old imports and plotting stuff.Here they are so you can count them:<jupyter_code>length, depth = 40, 100
model = 1 + np.tri(depth, length, -depth//3)
model[:depth//3,:] = 0
rocks = np.array([[2700, 2750], [2400, 2450], [2800, 3000]])
earth = np.take(rocks, model.astype(int), axis=0)
imp = np.apply_along_axis(np.product, -1, earth)
rc =  (imp[1:,:] - imp[:-1,:]) / (imp[1:,:] + imp[:-1,:])
w = bruges.filters.ricker(duration=0.100, dt=0.001, f=40)
synth = np.apply_along_axis(lambda t: np.convolve(t, w, mode='same'), axis=0, arr=rc)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11809171846128003, 'fraction_numerical': 0.02838904753905594, 'mean_word_length': 4.543761638733706, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 54, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '61832', 'license': 'Apache-2.0', 'n_tokens_mistral': 2086, 'n_tokens_neox': 1908, 'n_words': 777, 'path': 'notebooks/00_Synthetic_wedge_model.ipynb', 'repo_name': 'agile-geoscience/xlines', 'size': '5953'}