Looking at ways to build faulted models. We'll start with the prelims.
import matplotlib.pyplot as plt
import numpy as np
Set up the model parameters.
length = 100.# x range
width = 100. # y range
depth = 240. # z range
Define the plotting function.
def plot_slices(data, interpolation='nearest', cmap='Greys', vmin=0., vmax=1.):
plt.figure(figsize=(16,12))
for j in range(3):
plt.subplot(3,4,j+1)
plt.title('y = ' + str(int(j*(length/3))))
plt.imshow(data[:,:,int(j*(length/3))], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.subplot(3,4,4)
plt.title('y = ' + str(int(length-1)))
plt.imshow(data[:,:,-1], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
for j in range(3):
plt.subplot(3,4,5+j)
plt.title('x = ' + str(int(j*(width/3))))
plt.imshow(data[:,int(j*(width/3)),:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.subplot(3,4,8)
plt.title('x = ' + str(int(width-1)))
plt.imshow(data[:,-1,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
for j in range(3):
plt.subplot(3,4,9+j)
plt.title('z = ' + str(1 + int(1 + depth/3 + j*depth/(3*3))))
plt.imshow(data[1 + int(depth/3 + j*depth/(3*3)),:,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.subplot(3,4,12)
plt.title('z = ' + str(int(2*depth/3 - 1)))
plt.imshow(data[int(2*depth/3 - 1),:,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.show()
Let's start with thicknesses. We'll just normalize to a sum of 100% and split the total thickness accordingly. When assigning rocks, I guess we just cycle through them, bottom to top. Let's go straight to 3D.
# These are oldest > youngest
layers = [4,1,1,3,1,2,4]
# Start with 4 rocks
rocks = [(2800, 2850), (2700, 2750), (2400, 2450), (2300,2350)]
def make_ai(rocks):
rx = np.array(rocks)
return rx[:,0] * rx[:,1] / 10e6
ai = make_ai(rocks)
print ai
[ 0.798 0.7425 0.588 0.5405]
def get_depths(depth, layers):
l = [depth]
for layer in layers:
thickness = depth * float(layer) / sum(layers)
l.append(l[-1] - thickness)
return l[:-1]
depths = get_depths(depth, layers)
for i,d in enumerate(depths):
print ai[i%len(ai)]
0.798 0.7425 0.588 0.5405 0.798 0.7425 0.588
def make_earth(depth, length, width, layers, ai, throw=2):
'''
Given dimensions, a list of AIs, and a list of layers,
the function returns an array of AIs - a model earth.
'''
earth = np.empty((depth,length,width))
# Make a new set of layers with real thicknesses
# Actually want these as list of bases, bottom first
depths = get_depths(depth, layers)
for i,d in enumerate(depths):
earth[:d] = ai[i%len(ai)]
# Make the faulted layers (shifted version)
earth_f = np.empty_like(earth)
faulted_layers = layers[:]
faulted_layers[0] += throw
faulted_layers[-1] -= throw
depths = get_depths(depth, faulted_layers)
for i,d in enumerate(depths):
earth_f[:d] = ai[i%len(ai)]
earth[:,length/2:,:] = earth_f[:,length/2:,:]
return earth
Maybe there's a function to make a 'shifted' version of a vector?
%timeit make_earth(depth, length, width, layers, ai)
10 loops, best of 3: 34.6 ms per loop
earth = make_earth(depth, length, width, layers, ai)
plot_slices(earth, interpolation='nearest', cmap='jet')