X lines of Python

Wedge model

This is part of an Agile blog series called x lines of Python.

We start with the usual preliminaries.

Make an earth model

We'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 × n matrix with ones below the kth 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.

We'll have a quick look with some very basic plotting commands.

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 ith element from this other array, rocks in this case.

Edit: 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.

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.

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 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.

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.

We'll use a wavelet function from 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.

Let's make sure it looks OK:

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).

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: