This Notebook accompanies an Agile Geoscience blog post on the same subject.
A quick look at a progression of ways to parameterize a seismic attribute. We'll use the Canny filter, since I was playing with it anyway (see Sobel_filtering_horizons.ipynb).
Some ways to parameterize, in roughly increasing order of sophistication:
We could use all of the above approaches in a Notebook too — but why would we? There are better ways.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Next we need some data.
We will use the open-access Penobscot dataset from offshore Nova Scotia, quite near my house. This is a great dataset for lots of reasons, but especially because it is freely available and can be shared in almost any way you like. You can download it from the Open Seismic Repository. The data is owned by the Nova Scotia Deptartment of Energy and licensed CC-BY.
I saved this horizon as a NumPy object in the Notebook Sobel_filtering_horizons.ipynb.
horizon = np.load('Penobscot_HorB.npy')
plt.figure()
plt.imshow(-1 * horizon, aspect=0.5)
plt.show()
vmin = np.percentile(horizon[horizon > 0], 1)
vmax = np.percentile(horizon[horizon > 0], 99)
print 'min {0}, max {1}'.format(vmin, vmax)
We will compare our new attribute with this one later.
import scipy.ndimage
def do_sobel(image):
dx = scipy.ndimage.sobel(image, 0) # horizontal derivative
dy = scipy.ndimage.sobel(image, 1) # vertical derivative
mag = np.hypot(dx, dy) # magnitude
mag *= 255.0 / np.max(mag) # normalize
return mag
mag = do_sobel(horizon)
plt.imshow(mag, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=3)
plt.show()
plt.figure(figsize=(15,10))
plt.imshow(-1*horizon, aspect=0.5, cmap="cubehelix", vmin=-vmax, vmax=-vmin)
plt.colorbar(shrink=0.75) # shrink makes the colourbar a bit shorter
plt.show()
The nicest way to figure out how it works is to look at the example from SciKits website. Copy and paste and it runs perfectly...
# You might need to install this first: pip install scikit-image
from skimage import filter
# Generate noisy image of a square
im = np.zeros((128, 128))
im[32:-32, 32:-32] = 1
im = scipy.ndimage.rotate(im, 15, mode='constant')
im = scipy.ndimage.gaussian_filter(im, 4)
im += 0.2 * np.random.random(im.shape)
# Compute the Canny filter for two values of sigma
edges1 = filter.canny(im)
edges2 = filter.canny(im, sigma=3)
# display results
fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(12, 5))
ax1.imshow(im, cmap=plt.cm.jet)
ax1.axis('off')
ax1.set_title('noisy image', fontsize=20)
ax2.imshow(edges1, cmap=plt.cm.gray)
ax2.axis('off')
ax2.set_title('Canny filter, $\sigma=1$', fontsize=20)
ax3.imshow(edges2, cmap=plt.cm.gray)
ax3.axis('off')
ax3.set_title('Canny filter, $\sigma=3$', fontsize=20)
fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
bottom=0.02, left=0.02, right=0.98)
plt.show()
Now we can get to work. Let's apply some different ways to investigate the parameters of skimage.filter.canny
.
edge = filter.canny(horizon)
plt.figure(figsize=(15, 5))
plt.subplot(121)
plt.imshow(horizon, aspect=0.5, cmap='cubehelix_r', vmin=vmin, vmax=vmax)
plt.subplot(122)
plt.imshow(edge, aspect=0.5, cmap='Greys', vmin=0, vmax=1)
plt.show()