<jupyter_start><jupyter_text>Using a domain map to model a faultWe discussed the use of a binary map to delineate regions in which different potential fields are "active" - e.g. to model unconformities or faults. Let's have a look at a simple fault example here.<jupyter_code>import numpy as np
from scipy.interpolate import Rbf

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm

# for least square solution to line
from numpy.linalg import lstsq

plt.rcParams['figure.figsize'] = (8., 6.)
plt.rcParams['font.size'] = 16.<jupyter_output><empty_output><jupyter_text>Create randomized dataAs before, we will include a data generation procedure, in order to be able to test the effect of varying resolution on model results below.<jupyter_code>n_pts = 10 # Points per layer

# Linear functions for line data
l1 = lambda x : 0.25*x + 30
l2 = lambda x : 0.25*x + 40
l3 = lambda x : 0.25*x + 50

# set seed for reproducibility
np.random.seed(123)
# sampling points
l1_pts_x = np.random.uniform(0,90,n_pts)
l1_pts_y = l1(l1_pts_x)
l2_pts_x = np.random.uniform(0,90,n_pts)
l2_pts_y = l2(l2_pts_x)
l3_pts_x = np.random.uniform(0,90,n_pts)
l3_pts_y = l3(l3_pts_x)<jupyter_output><empty_output><jupyter_text>Define faultIn order to define a fault, we will model it here as a line through two points:<jupyter_code>fault_point_1 = (40,60)
fault_point_2 = (60,20)
x_coords, y_coords = [40,60], [60,20] # zip(*points)
A = np.vstack([x_coords, np.ones(len(x_coords))]).T
m, c = lstsq(A, y_coords, rcond=None)[0]
print("Line Solution is y = {m}x + {c}".format(m=m,c=c))<jupyter_output>Line Solution is y = -2.0000000000000013x + 140.00000000000014<jupyter_text>We next define the effect of this fault on the surrounding points as an offset:<jupyter_code>offset = 10 # offset of block on right side of fault
f = lambda x : m*x + c

# Create filters to determine points on each side of fault
filter_l1 = f(l1_pts_x) < l1_pts_y
filter_l2 = f(l2_pts_x) < l2_pts_y
filter_l3 = f(l3_pts_x) < l3_pts_y

# create copies of arrays to avoid confusion...
l1_pts_x_fault = l1_pts_x.copy()
l1_pts_y_fault = l1_pts_y.copy()
l2_pts_x_fault = l2_pts_x.copy()
l2_pts_y_fault = l2_pts_y.copy()
l3_pts_x_fault = l3_pts_x.copy()
l3_pts_y_fault = l3_pts_y.copy()

# Adjust y-values
l1_pts_y_fault[filter_l1] -= offset
l2_pts_y_fault[filter_l2] -= offset
l3_pts_y_fault[filter_l3] -= offset
# Adjust x-values
l1_pts_x_fault[filter_l1] -= 1/m*offset
l2_pts_x_fault[filter_l2] -= 1/m*offset
l3_pts_x_fault[filter_l3] -= 1/m*offset
f = lambda x : m*x + c

plt.plot(l1_pts_x_fault, l1_pts_y_fault, 'o')
plt.plot(l2_pts_x_fault, l2_pts_y_fault, 'o')
plt.plot(l3_pts_x_fault, l3_pts_y_fault, 'o')

plt.plot(x_coords, y_coords, 'ro')
xrange = np.arange(100)
plt.plot(xrange, f(xrange))

plt.ylim([0,100])

# plt.axis('equal')<jupyter_output><empty_output><jupyter_text>_Note_: this is actually a very primitive definition of a line: we use here a so-called "kinematic model" as we explicitly define the fault position and the effect on the surrounding data.  Create domain map for fault Interpolate scalar field for fault domainsAs a first step, we will create a domain map on the basis of the fault. Towards this aim, we will also generate a scalar field interpolation for the fault itself:<jupyter_code># set up grid
ti = np.linspace(0, 100, 100)
XI, YI = np.meshgrid(ti, ti)
# trick: add points at lower left and upper right to obtain gradients in field
# (otherwise, we would simply interpolate zeros...)
x_coords_ext = np.hstack([x_coords, 0, 100])
y_coords_ext = np.hstack([y_coords, 0, 100])
rbf = Rbf(x_coords_ext, y_coords_ext, [0,0,-1,1], function="gaussian", epsilon=1000)
ZI = rbf(XI, YI)
plt.contour(ZI, 0)
plt.plot(x_coords, y_coords, 'ro')
# plt.colorbar()<jupyter_output><empty_output><jupyter_text>Extract domain mapIn this example, we have two domains, one on each side of the fault. In this case, we have a simple boolean map (for generalizations, we can extend it to an integer map for multiple domains):<jupyter_code>fault_domains = np.zeros_like(ZI)
fault_domains[ZI>0] = 1
plt.pcolor(fault_domains, cmap='RdBu', alpha=0.2)
plt.contour(ZI, 0)
plt.plot(x_coords, y_coords, 'ro')<jupyter_output><empty_output><jupyter_text>Interpolate fields on both side of faultNow that we have the domains defined, we can interpolate the scalar fields on both sides of the fault. We will here use a simple "divide-and-conquer" method, i.e.:- separate point data sets according to domain in which points are in (note: we can also very effectively use the fault scalar field for this purpose);- interpolate a separate scalar field for each domain;- finally, create one combined scalar field with values according to each domain;- extract isosurfaces using contouring method (marching squares). Separate points on both sidesLet's start with the separation of points:<jupyter_code># left domain
l1_pts_x_left = l1_pts_x_fault[np.where(rbf(l1_pts_x_fault, l1_pts_y_fault) < 0)]
l1_pts_y_left = l1_pts_y_fault[np.where(rbf(l1_pts_x_fault, l1_pts_y_fault) < 0)]
l2_pts_x_left = l2_pts_x_fault[np.where(rbf(l2_pts_x_fault, l2_pts_y_fault) < 0)]
l2_pts_y_left = l2_pts_y_fault[np.where(rbf(l2_pts_x_fault, l2_pts_y_fault) < 0)]
l3_pts_x_left = l3_pts_x_fault[np.where(rbf(l3_pts_x_fault, l3_pts_y_fault) < 0)]
l3_pts_y_left = l3_pts_y_fault[np.where(rbf(l3_pts_x_fault, l3_pts_y_fault) < 0)]
# right domain
l1_pts_x_right = l1_pts_x_fault[np.where(rbf(l1_pts_x_fault, l1_pts_y_fault) > 0)]
l1_pts_y_right = l1_pts_y_fault[np.where(rbf(l1_pts_x_fault, l1_pts_y_fault) > 0)]
l2_pts_x_right = l2_pts_x_fault[np.where(rbf(l2_pts_x_fault, l2_pts_y_fault) > 0)]
l2_pts_y_right = l2_pts_y_fault[np.where(rbf(l2_pts_x_fault, l2_pts_y_fault) > 0)]
l3_pts_x_right = l3_pts_x_fault[np.where(rbf(l3_pts_x_fault, l3_pts_y_fault) > 0)]
l3_pts_y_right = l3_pts_y_fault[np.where(rbf(l3_pts_x_fault, l3_pts_y_fault) > 0)]


fit, axes = plt.subplots(1,2, figsize=(14,6))

axes[0].pcolor(fault_domains, cmap='RdBu', alpha=0.2)
axes[0].contour(ZI, 0)
axes[0].plot(x_coords, y_coords, 'ro')
# plot points on left side
axes[0].plot(l1_pts_x_left, l1_pts_y_left, 'o')
axes[0].plot(l2_pts_x_left, l2_pts_y_left, 'o')
axes[0].plot(l3_pts_x_left, l3_pts_y_left, 'o')
axes[0].set_title("Points in left domain")

axes[1].pcolor(fault_domains, cmap='RdBu', alpha=0.2)
axes[1].contour(ZI, 0)
axes[1].plot(x_coords, y_coords, 'ro')
# plot points on right side
axes[1].set_title("Points in right domain")
axes[1].plot(l1_pts_x_right, l1_pts_y_right, 'o')
axes[1].plot(l2_pts_x_right, l2_pts_y_right, 'o')
axes[1].plot(l3_pts_x_right, l3_pts_y_right, 'o');<jupyter_output><empty_output><jupyter_text>Interpolate scalar fields for each set of pointsWe now interpolate one scalar field for the set of points in each domain:<jupyter_code># left domain

# combine data in arrays
x = np.hstack([l1_pts_x_left, l2_pts_x_left, l3_pts_x_left])
y = np.hstack([l1_pts_y_left, l2_pts_y_left, l3_pts_y_left])

# give points values
z = np.hstack([np.ones(len(l1_pts_x_left))*10, 
               np.ones(len(l2_pts_x_left))*20, 
               np.ones(len(l3_pts_x_left))*30])

# use RBF
rbf = Rbf(x, y, z, function="multiquadric", epsilon=50)
ZI_left = rbf(XI, YI)
# right domain

# combine data in arrays
x = np.hstack([l1_pts_x_right, l2_pts_x_right, l3_pts_x_right])
y = np.hstack([l1_pts_y_right, l2_pts_y_right, l3_pts_y_right])

# give points values
z = np.hstack([np.ones(len(l1_pts_x_right))*10, 
               np.ones(len(l2_pts_x_right))*20, 
               np.ones(len(l3_pts_x_right))*30])

# use RBF
rbf = Rbf(x, y, z, function="multiquadric", epsilon=50)
ZI_right = rbf(XI, YI)<jupyter_output><empty_output><jupyter_text>Let's create a plot for both interpolated fields:<jupyter_code>fit, axes = plt.subplots(1,3, figsize=(16,5))

# Left field

axes[0].pcolor(ZI_left, cmap='viridis', alpha=0.2)
axes[0].contour(ZI_left, [10,20,30])
axes[0].contour(ZI, 0)
axes[0].plot(x_coords, y_coords, 'ro')
# plot points on left side
axes[0].plot(l1_pts_x_left, l1_pts_y_left, 'o')
axes[0].plot(l2_pts_x_left, l2_pts_y_left, 'o')
axes[0].plot(l3_pts_x_left, l3_pts_y_left, 'o')
axes[0].set_title("Points in left domain")

# Right field

axes[1].pcolor(ZI_right, cmap='viridis', alpha=0.2)
axes[1].contour(ZI_right, [10,20,30])
axes[1].contour(ZI, 0)
axes[1].plot(x_coords, y_coords, 'ro')
# plot points on right side
axes[1].set_title("Points in right domain")
axes[1].plot(l1_pts_x_right, l1_pts_y_right, 'o')
axes[1].plot(l2_pts_x_right, l2_pts_y_right, 'o')
axes[1].plot(l3_pts_x_right, l3_pts_y_right, 'o');

# Combined field

axes[2].pcolor(ZI_comb, cmap='viridis', alpha=0.2)
axes[2].contour(ZI_comb, [10,20,30])
axes[2].contour(ZI, 0)
axes[2].plot(x_coords, y_coords, 'ro')
# plot all points
axes[2].set_title("Combined field and all points")
axes[2].plot(l1_pts_x_fault, l1_pts_y_fault, 'o')
axes[2].plot(l2_pts_x_fault, l2_pts_y_fault, 'o')
axes[2].plot(l3_pts_x_fault, l3_pts_y_fault, 'o')

fit, axes = plt.subplots(1,2, figsize=(14,6))

axes[0].pcolor(ZI_left, cmap='viridis', alpha=0.2)
axes[0].contour(ZI_left, [10,20,30])
axes[0].contour(ZI, 0)
axes[0].plot(x_coords, y_coords, 'ro')
# plot points on left side
axes[0].plot(l1_pts_x_left, l1_pts_y_left, 'o')
axes[0].plot(l2_pts_x_left, l2_pts_y_left, 'o')
axes[0].plot(l3_pts_x_left, l3_pts_y_left, 'o')
axes[0].set_title("Points in left domain")

axes[1].pcolor(ZI_right, cmap='viridis', alpha=0.2)
axes[1].contour(ZI_right, [10,20,30])
axes[1].contour(ZI, 0)
axes[1].plot(x_coords, y_coords, 'ro')
# plot points on right side
axes[1].set_title("Points in right domain")
axes[1].plot(l1_pts_x_right, l1_pts_y_right, 'o')
axes[1].plot(l2_pts_x_right, l2_pts_y_right, 'o')
axes[1].plot(l3_pts_x_right, l3_pts_y_right, 'o');<jupyter_output><empty_output><jupyter_text>Combine scalar fieldsIf we chose a reasonable interpolation, then we now have interpolated scalar fields with reasonable fits on both sides of the fault. In the next step, we will combine both fields, according to the (fault) domain map. Note that this simply means that we use map the scalar field values according to each domain:<jupyter_code>ZI_comb = np.zeros_like(ZI_left)

ZI_comb[ZI < 0] = ZI_left[ZI < 0]
ZI_comb[ZI > 0] = ZI_right[ZI > 0]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12828272966904228, 'fraction_numerical': 0.043053792834130626, 'mean_word_length': 6.33810888252149, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '84042', 'license': 'MIT', 'n_tokens_mistral': 4667, 'n_tokens_neox': 4339, 'n_words': 977, 'path': 'geomodeling/domain_maps_fault.ipynb', 'repo_name': 'cgre-aachen/teaching', 'size': '10243'}