<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Simple model
This notebook is a series of independent cells showing how to create a simple model from the beginning to the end using GeMpy
Step1:  Importing dependencies
Step3:  Visualize data
Step4:  Interactive pandas Dataframe
Using qgrid it is possible to modify the tables in place as following
Step5:  Grid and potential field
We can see the potential field generated out of the data above
Step6:  From potential field to block
The potential field describe the deposition form and direction of a basin. However, in most scenarios the real goal of structural modeling is the segmentation in layers of areas with significant change of properties (e.g. shales and carbonates). Since we need to provide at least one point per interface, we can easily compute the value of the potential field at the intersections between two layers. Therefore, by simple comparison between a concrete value of the potential field and the values of the interfaces it is possible to segment the domain into layers Fig X.
Step8:  Combining potential fields
Step9:  This potential field gives the following block
Step10:  Combining both potential field where the first potential field is younger than the second we can obtain the following structure.
Step11:  Side note
<ASSISTANT_TASK:>
Python Code:
from IPython.display import HTML
HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Explanation: Simple model
This notebook is a series of independent cells showing how to create a simple model from the beginning to the end using GeMpy
End of explanation
# Importing
import importlib
import theano.tensor as T
import sys, os
sys.path.append("/home/bl3/PycharmProjects/GeMpy/GeMpy")
import pandas as pn
import GeMpy_core
import Visualization
import numpy as np
importlib.reload(GeMpy_core)
os.environ['CUDA_LAUNCH_BLOCKING'] = '1'
np.set_printoptions(precision = 4, linewidth= 300, suppress =  True)
import matplotlib.pyplot as plt
from matplotlib import cm
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'png')
plt.rcParams['savefig.dpi'] = 75
plt.rcParams['figure.autolayout'] = False
plt.rcParams['figure.figsize'] = 10, 6
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['axes.titlesize'] = 20
plt.rcParams['font.size'] = 16
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.markersize'] = 8
plt.rcParams['legend.fontsize'] = 14
import seaborn as sns
# This sets reasonable defaults for font size for
# a figure that will go in a paper
sns.set_context("paper")
# Set the font to be serif, rather than sans
sns.set(font='Arial')
# Make the background white, and specify the
# specific font family
sns.set_style("white")
import qgrid
qgrid.nbinstall(overwrite=True)  # copies javascript dependencies to your /nbextensions folder
%matplotlib inline
Explanation: Importing dependencies
End of explanation
# Setting extend, grid and compile
# Setting the extent
test = GeMpy_core.GeMpy()
test.import_data([0,10,0,10,0,10])
# =========================
# DATA GENERATION IN PYTHON
# =========================
# Layers coordinates
layer_1 = np.array([[0.5,4,7], [2,4,6.5], [4,4,7], [5,4,6]])#-np.array([5,5,4]))/8+0.5
layer_2 = np.array([[3,4,5], [6,4,4],[8,4,4], [7,4,3], [1,4,6]])
layers = np.asarray([layer_1,layer_2])
# Foliations coordinates
dip_pos_1 = np.array([7,4,7])#- np.array([5,5,4]))/8+0.5
dip_pos_2 = np.array([2.,4,4])
# Dips
dip_angle_1 = float(15)
dip_angle_2 = float(340)
dips_angles = np.asarray([dip_angle_1, dip_angle_2], dtype="float64")
# Azimuths
azimuths = np.asarray([90,90], dtype="float64")
# Polarity
polarity = np.asarray([1,1], dtype="float64")
# Pandas Dataframe with the interfaces data
test.Data.Interfaces = pn.DataFrame(
    data = {"X" :np.append(layer_1[:, 0],layer_2[:,0]),
            "Y" :np.append(layer_1[:, 1],layer_2[:,1]),
            "Z" :np.append(layer_1[:, 2],layer_2[:,2]),
            "formation" : np.append(
               np.tile("Layer 1", len(layer_1)), 
               np.tile("Layer 2", len(layer_2)))})
# Pandas Dataframe with the Foliations data
test.Data.Foliations =   pn.DataFrame(
    data = {"X" :np.append(dip_pos_1[0],dip_pos_2[0]),
            "Y" :np.append(dip_pos_1[ 1],dip_pos_2[1]),
            "Z" :np.append(dip_pos_1[ 2],dip_pos_2[2]),
            "azimuth" : azimuths,
            "dip" : dips_angles,
            "polarity" : polarity,
            "formation" : ["Layer 1", "Layer 2"]})
# Creation of the formations (to be deprecated)
test.Data.formations = test.Data.Interfaces["formation"].unique()
# Calculation of the pola gradient from dip azimuth and polarity
test.Data.calculate_gradient()
# Set defautl series
test.Data.set_series()
# Method to be sure all objects get updated
test.update_data()
# The following code is for the visualization of labels of the input data (yet to be implemented)
# ----------------------------------------------------------------------------------------------
def annotate_plot(frame, label_col, x, y, **kwargs):
    
    Annotate the plot of a given DataFrame using one of its columns
    Should be called right after a DataFrame or series plot method,
    before telling matplotlib to show the plot.
    Parameters
    ----------
    frame : pandas.DataFrame
    plot_col : str
        The string identifying the column of frame that was plotted
    label_col : str
        The string identifying the column of frame to be used as label
    kwargs:
        Other key-word args that should be passed to plt.annotate
    Returns
    -------
    None
    Notes
    -----
    After calling this function you should call plt.show() to get the
    results. This function only adds the annotations, it doesn't show
    them.
    
    import matplotlib.pyplot as plt  # Make sure we have pyplot as plt
    for label, x, y in zip(frame[label_col], frame[x], frame[y]):
        plt.annotate(label, xy=(x+0.2, y+0.15), **kwargs)
        
inter_labels =[r'${\bf{x}}_{\alpha \, 0}^1$',
               r'${\bf{x}}_{\alpha \, 1}^1$',
               r'${\bf{x}}_{\alpha \, 2}^1$',
               r'${\bf{x}}_{\alpha \, 3}^1$',
               r'${\bf{x}}_{\alpha \, 0}^2$',
               r'${\bf{x}}_{\alpha \, 1}^2$',
               r'${\bf{x}}_{\alpha \, 2}^2$',
               r'${\bf{x}}_{\alpha \, 3}^2$',
               r'${\bf{x}}_{\alpha \, 4}^2$']              
        
foli_labels =[r'${\bf{x}}_{\beta \,{0}}$',
              r'${\bf{x}}_{\beta \,{1}}$']        
test.Data.Interfaces['labels'] = pn.Series(inter_labels)
test.Data.Foliations['labels'] = pn.Series(foli_labels)
# Plot and table
test.Plot.plot_data()
annotate_plot(test.Data.Interfaces, 'labels','X', 'Z', size = 'x-large')
annotate_plot(test.Data.Foliations, 'labels','X', 'Z', size = 'x-large')
test.Data.Interfaces
Explanation: Visualize data
End of explanation
import qgrid
qgrid.show_grid(test.Data.Interfaces)
Explanation: Interactive pandas Dataframe
Using qgrid it is possible to modify the tables in place as following:
End of explanation
test.create_grid()
test.set_interpolator(u_grade=0, verbose = 0)
test.Plot.plot_potential_field(4, direction='y', colorbar = True, cmap = 'magma')
#test.Interpolator.DK;
from ipywidgets import widgets
from ipywidgets import interact
def cov_cubic_f(r,a = 6, c_o = 1):
    if r <= a:
        return c_o*(1-7*(r/a)**2+35/4*(r/a)**3-7/2*(r/a)**5+3/4*(r/a)**7)
    else:
        return 0
    
def cov_cubic_d1_f(r,a = 6., c_o = 1):
    SED_dips_dips = r
    f = c_o
    return  (f * ((-14 /a ** 2) + 105 / 4 * SED_dips_dips / a ** 3 -
            35 / 2 * SED_dips_dips ** 3 / a ** 5 + 21 / 4 * SED_dips_dips ** 5 / a ** 7))
def cov_cubic_d2_f(r, a = 6, c_o = 1):
    SED_dips_dips = r
    f = c_o
    return 7*f*(9*r**5-20*a**2*r**3+15*a**4*r-4*a**5)/(2*a**7)
def plot_potential_var(a = 10, c_o = 1, nugget_effect = 0):
    x = np.linspace(0,12,50)
    y = [cov_cubic_f(i, a = a, c_o = c_o) for i in x]
    
    fig = plt.figure()
    ax1 = fig.add_subplot(121)
    ax1.plot(x,c_o-np.asarray(y)+nugget_effect)
    plt.hlines(0,0,12, linestyles = "--")
    plt.title("Variogram")
    plt.margins(0,0.1)
    
    ax2 = fig.add_subplot(122)
    ax2.plot(x,np.asarray(y))
    ax2.scatter(0,nugget_effect+c_o)
    plt.title("Covariance Function")
    plt.tight_layout()
    plt.margins(0,0.1)
    plt.suptitle('$C_Z(r)$', y = 1.08, fontsize=15, fontweight='bold')
def plot_potential_direction_var( a = 10, c_o = 1, nugget_effect = 0):
    x = np.linspace(0,12,50)
    y = np.asarray([cov_cubic_d1_f(i, a = a, c_o = c_o) for i in x])
    
    fig = plt.figure()
    ax1 = fig.add_subplot(121)
    ax1.plot(x,c_o-np.asarray(y)+nugget_effect)
    plt.title("Variogram")
    plt.margins(0,0.1)
    
    ax2 = fig.add_subplot(122)
    ax2.plot(x,np.asarray(y))
    #ax2.scatter(0,c_o)
    plt.title("Cross-Covariance Function")
    plt.tight_layout()
    plt.margins(0,0.1)
    plt.suptitle('$C\'_Z / r$', y = 1.08, fontsize=15, fontweight='bold')    
def plot_directionU_directionU_var(a = 10, c_o = 1, nugget_effect = 0):
    
    x = np.linspace(0.01,12,50)
    d1 = np.asarray([cov_cubic_d1_f(i, a = a, c_o = c_o) for i in x])
    d2 = np.asarray([cov_cubic_d2_f(i, a = a, c_o = c_o) for i in x])
    
    y = -(d2) # (0.5*x**2)/(x**2)*
    
    fig = plt.figure()
    ax1 = fig.add_subplot(121)
    ax1.plot(x,c_o-np.asarray(y)+nugget_effect)
    plt.title("Variogram")
    plt.margins(0,0.1)
    
    ax2 = fig.add_subplot(122)
    ax2.plot(x,np.asarray(y))
    ax2.scatter(0,nugget_effect+y[0], s = 20)
    plt.title("Covariance Function")
    plt.tight_layout()
    plt.margins(0,0.1)
    plt.suptitle('$C_{\partial {Z}/ \partial x, \, \partial {Z}/ \partial x}(h_x)$'
                 , y = 1.08, fontsize=15)
    
def plot_directionU_directionV_var(a = 10, c_o = 1, nugget_effect = 0):
    x = np.linspace(0.01,12,50)
    d1 = np.asarray([cov_cubic_d1_f(i, a = a, c_o = c_o) for i in x])
    d2 = np.asarray([cov_cubic_d2_f(i, a = a, c_o = c_o) for i in x])
    y = -(d2-d1) # (0.5*x**2)/(x**2)*
    
    fig = plt.figure()
    ax1 = fig.add_subplot(121)
    ax1.plot(x,c_o-np.asarray(y)+nugget_effect)
    plt.title("Variogram")
    plt.margins(0,0.1)
    
    ax2 = fig.add_subplot(122)
    ax2.plot(x,np.asarray(y))
    ax2.scatter(0,nugget_effect+y[0], s = 20)
    plt.title("Covariance Function")
    plt.tight_layout()
    plt.margins(0,0.1)
    plt.suptitle('$C_{\partial {Z}/ \partial x, \, \partial {Z}/ \partial y}(h_x,h_y)$'
                 , y = 1.08, fontsize=15)
    
def plot_all(a = 10, c_o = 1, nugget_effect = 0):
    
    
    plot_potential_direction_var(a, c_o, nugget_effect)
    plot_directionU_directionU_var(a, c_o, nugget_effect)
    plot_directionU_directionV_var(a, c_o, nugget_effect)
Explanation: Grid and potential field
We can see the potential field generated out of the data above
End of explanation
test.Interpolator.compute_block_model()
test.Plot.plot_block_section()
Explanation: From potential field to block
The potential field describe the deposition form and direction of a basin. However, in most scenarios the real goal of structural modeling is the segmentation in layers of areas with significant change of properties (e.g. shales and carbonates). Since we need to provide at least one point per interface, we can easily compute the value of the potential field at the intersections between two layers. Therefore, by simple comparison between a concrete value of the potential field and the values of the interfaces it is possible to segment the domain into layers Fig X.
End of explanation
# Setting extend, grid and compile
# Setting the extent
two_pot = GeMpy_core.GeMpy()
two_pot.import_data([0,10,0,10,0,10])
# Data Generation
layer_1 = np.array([[0.5,4,7], [2,4,6.5], [4,4,7], [5,4,6]])#-np.array([5,5,4]))/8+0.5
layer_2 = np.array([[3,4,5], [6,4,4],[8,4,4], [7,4,3], [1,4,6]])
layer_3 = np.array([[2,4,3], [8,4,2], [9,4,3]])
dip_pos_1 = np.array([7,4,7])#- np.array([5,5,4]))/8+0.5
dip_pos_2 = np.array([2.,4,4])
dip_pos_3 = np.array([1,4,1])
dip_angle_1 = float(15)
dip_angle_2 = float(340)
dip_angle_3 = float(80)
dips_angles = np.asarray([dip_angle_1, dip_angle_2, dip_angle_3], dtype="float64")
layers = np.asarray([layer_1,layer_2, layer_3])
azimuths = np.asarray([90,90, 90], dtype="float64")
polarity = np.asarray([1,1, 1], dtype="float64")
two_pot.Data.Interfaces = pn.DataFrame(data = {"X" :np.hstack((layer_1[:, 0],layer_2[:,0], layer_3[:,0])),
                                               "Y" :np.hstack((layer_1[:, 1],layer_2[:,1], layer_3[:,1])),
                                               "Z" :np.hstack((layer_1[:, 2],layer_2[:,2], layer_3[:,2])),
"formation" : np.hstack((np.tile("Layer 1", len(layer_1)), np.tile("Layer 2", len(layer_2)),
                        np.tile("Layer 3", len(layer_3))))})
two_pot.Data.Foliations =   pn.DataFrame(data = {"X" :np.hstack((dip_pos_1[0], dip_pos_2[0], dip_pos_3[0])),
                                                 "Y" :np.hstack((dip_pos_1[ 1],dip_pos_2[1], dip_pos_3[1])),
                                                 "Z" :np.hstack((dip_pos_1[ 2],dip_pos_2[2], dip_pos_3[2])),
                                         "azimuth" : azimuths,
                                         "dip" : dips_angles,
                                         "polarity" : polarity,
"formation" : ["Layer 1", "Layer 2", 'Layer 3']})
inter_labels =[r'${\bf{x}}_{\alpha \, 0}^1$',
               r'${\bf{x}}_{\alpha \, 1}^1$',
               r'${\bf{x}}_{\alpha \, 2}^1$',
               r'${\bf{x}}_{\alpha \, 3}^1$',
               r'${\bf{x}}_{\alpha \, 0}^2$',
               r'${\bf{x}}_{\alpha \, 1}^2$',
               r'${\bf{x}}_{\alpha \, 2}^2$',
               r'${\bf{x}}_{\alpha \, 3}^2$',
               r'${\bf{x}}_{\alpha \, 4}^2$',
               r'${\bf{x}}_{\alpha \, 0}^1$',
               r'${\bf{x}}_{\alpha \, 1}^1$',
               r'${\bf{x}}_{\alpha \, 2}^1$']              
        
foli_labels =[r'${\bf{x}}_{\beta \,{0}}$',
              r'${\bf{x}}_{\beta \,{1}}$',
              r'${\bf{x}}_{\beta \,{0}}$']
two_pot.Data.Interfaces['labels'] = pn.Series(inter_labels)
two_pot.Data.Foliations['labels'] = pn.Series(foli_labels)
two_pot.Data.formations = two_pot.Data.Interfaces["formation"].unique()
two_pot.Data.calculate_gradient()
two_pot.Data.set_series({'younger': ('Layer 1', 'Layer 2'),
                         'older': 'Layer 3'}, order = ['younger', 'older'])
two_pot.update_data()
def annotate_plot(frame, label_col, x, y, **kwargs):
    
    Annotate the plot of a given DataFrame using one of its columns
    Should be called right after a DataFrame or series plot method,
    before telling matplotlib to show the plot.
    Parameters
    ----------
    frame : pandas.DataFrame
    plot_col : str
        The string identifying the column of frame that was plotted
    label_col : str
        The string identifying the column of frame to be used as label
    kwargs:
        Other key-word args that should be passed to plt.annotate
    Returns
    -------
    None
    Notes
    -----
    After calling this function you should call plt.show() to get the
    results. This function only adds the annotations, it doesn't show
    them.
    
    import matplotlib.pyplot as plt  # Make sure we have pyplot as plt
    for label, x, y in zip(frame[label_col], frame[x], frame[y]):
        plt.annotate(label, xy=(x+0.2, y+0.15), **kwargs)
        
serie_to_plot = 'older'        
two_pot.Plot.plot_data(series = serie_to_plot)
annotate_plot(two_pot.Data.Interfaces[two_pot.Data.Interfaces['series'] == serie_to_plot]
              , 'labels','X', 'Z', size = 'x-large')
annotate_plot(two_pot.Data.Foliations[two_pot.Data.Foliations['series'] == serie_to_plot],
              'labels','X', 'Z', size = 'x-large')
two_pot.create_grid()
two_pot.set_interpolator(u_grade=0, verbose = 0)
Explanation: Combining potential fields: Depositional series
In reality, most geological settings are formed by a concatenation of depositional phases separated clearly by unconformity bounderies. Each of these phases can be model by a potential field. In order to capture this behavior, we can classify the formations that belong to individual depositional phase into categories or series. The potential field computed for each of these series could be seen as a sort of evolution of the basin if an unconformity did not occur. Finally, sorting the temporal relation between series allow to superpose the corresponding potential field at an specific location.
In the next example, we add a new serie consisting in a layer---'Layer 3'--- Fig X, which generate the potential field of Fig X and subsequently the block Figure X.
End of explanation
two_pot.Plot.plot_potential_field(4, n_pf=1, direction='y', colorbar = True, cmap = 'magma')
# Reset the block
two_pot.Interpolator.block.set_value(np.zeros_like(two_pot.Grid.grid[:,0]))
two_pot.Interpolator.compute_block_model(series_number=[1])
two_pot.Plot.plot_block_section()
Explanation: This potential field gives the following block
End of explanation
# Reset the block
two_pot.Interpolator.block.set_value(np.zeros_like(two_pot.Grid.grid[:,0]))
two_pot.Interpolator.compute_block_model(series_number=[0,1])
two_pot.Plot.plot_block_section()
Explanation: Combining both potential field where the first potential field is younger than the second we can obtain the following structure.
End of explanation
plot_potential_var(10,10**2 / 14 / 3  , 0.01)
plot_all(10,10**2 / 14 / 3  , 0.01) # 0**2 /14/3
Explanation: Side note: Example of covariances involved in the cokriging system
End of explanation
<END_TASK>