<jupyter_start><jupyter_text>Building geological models with LoopStructuralThis notebook will demonstrate how to build a basic geological model using LoopStructural.You will:* Understand how to format and generate an input dataset for LoopStructural* Add faults to the model* Add unconformities to the model* Visualise the model* Export the modelAdditional resources:* [LoopStructural Documentation](https://loop3d.github.io/LoopStructural/)  Implicit functionThe aim of implicit modelling is to find a function which represents the distance away from a reference geological horizon.There is no known analytical solution for this problem which means that they need to be approximated from the observations that are provided.We can approximate the implicit function using a weighted combination of basis functions: $f(x,y,z) = \sum^N_{i=0} v_i \varphi_i(X) $Geological observations including the location of contacts and orientation of stratigraphy can be converted into mathematical expressions:* Observations constraining the value of the scalar field     * $f(x,y,z) = v$* Observations constraining an form surface or interface between stratigraphic units     * $ \sum^N_{i=0} \sum^N_{j=i} f(x_i,y_i,z_i) - f(x_j,y_j,z_j) = 0$* Observations constraining the magnitude of the gradient norm of the scalar field       * $ \nabla f(x,y,z) = \textbf{n}$* Observations constraining a vector to be orthogonal to the scalar field    * $\nabla f(x,y,z) \cdot \textbf{t} = 0$      Imports for LoopStructural* To build a geological model using LoopStructural this can be done using only the `GeologicalModel` class.* To visualise the geologial model, there is a wrapper around lavavu that allows for all of the objects in the GeologicalModel to be visualised.<jupyter_code>from LoopStructural import GeologicalModel
from LoopStructural.visualisation import LavaVuModelViewer

from LoopStructural.datasets import load_claudius 
data, bb = load_claudius()<jupyter_output><empty_output><jupyter_text>LoopStructural data formatting* X - x component of the cartesian coordinates* Y - y component of the cartesian coordinates* Z - z component of the cartesian coordinates* feature_name - unique name of the geological feature being modelled* val - value observations of the scalar field* interface - unique identifier for an inteface containing similar scalar field values* nx - x component of the gradient norm* ny - y component of the gradient norm* nz - z component of the gradient norm* gx - x component of a gradient constraint* gy - y component of a gradient constraint* gz - z component of a gradient constraint* tx - x component of a gradient tangent constraint* ty - y component of a gradient tangent constraint* tz - z component of a gradient tangent constraint* coord - coordinate of the structural frame data point is used for<jupyter_code>data<jupyter_output><empty_output><jupyter_text>GeologicalModel Define the model bounding box using an origin and maximum vector. The model will be reprojected so that the origin is 0,0,0. The model can be rescaled so that the maximum length is 1 using the `rescale` flag.<jupyter_code>model = GeologicalModel(bb[0,:],bb[1,:],rescale=False)<jupyter_output><empty_output><jupyter_text>Link the model to the pandas DataFrame, during this set the model data will be rescaled/reprojected.<jupyter_code>model.set_model_data(data)<jupyter_output><empty_output><jupyter_text>LoopStructural uses a time aware approach building the geological model. This means that features in the model are added backwards in time with the youngest feature being added first. To add features to the model there are different methods that are used depending on the type of feature being modelled.* `create_and_add_foliation` is used to add a foliation (e.g. stratigraphy)* `create_and_add_fault`If there are unconformities, the default LoopStructural behaviour is to add an unconformity at the $0$ isovalue for the younger feature.<jupyter_code>model.create_and_add_foliation('strati')<jupyter_output><empty_output><jupyter_text>LoopStructural assembles the linear systems and determines the relationship between geological features but does not solve the linear systems. To solve the geological model, the `update` method needs to be called.<jupyter_code>model.update()<jupyter_output><empty_output><jupyter_text>Each `GeologicalFeature` can be accessed from the `GeologicalModel` by either accessing the feature directly `model.features[0]`, or using the index operator with the feature name e.g. `model['strati']`.The `GeologicalFeature` can be evaluated for the value or gradient of the scalar field* `evaluate_value(xyz)` which will return the scalar field value for the locations xyz* `evaluate_gradient(xyz)` which will return the gradient of the scalar field at the locations xyzBoth will return `np.nan` where the scalar field is not defined.  Visualising a model`LavaVuModelViewer` is used for visualising a LoopStructural model, it is a wrapper for lavavu and automatically processes the geological model for visualisation. A guide for the visualisation using LoopStructural is [here](https://loop3d.github.io/LoopStructural/user_guide/visualisation.html)<jupyter_code>view = LavaVuModelViewer(model)
view.add_isosurface(model['strati'],slices=data['val'].unique())
view.add_data(model['strati'])
view.interactive()<jupyter_output>2021-04-20 01:56:34,114 ~ LoopStructural.visualisation.model_visualisation ~ WARNING    ~ Cannot isosurface strati at nan, skipping
2021-04-20 01:56:34,114 ~ LoopStructural.visualisation.model_visualisation ~ WARNING    ~ Cannot isosurface strati at nan, skipping
2021-04-20 01:56:34,114 ~ LoopStructural.visualisation.model_visualisation ~ WARNING    ~ Cannot isosurface strati at nan, skipping<jupyter_text>An example with faults and an unconformity  Generating a dataset for LoopStructuralWe will define two stratigraphic groups:1. upper group which is defined by two horizons that are undulating and is unconformably overlies the lower group 2. lower group a flat surface that is faulted by a reverse fault  Define the upper group by randomly generating two pointsets for a value of 0 and 1<jupyter_code>import numpy as np
xx,yy = np.meshgrid(np.linspace(0,10,100),np.linspace(0,10,100))
pts = np.vstack([xx.flatten(),yy.flatten(),np.zeros_like(xx.flatten()),np.zeros_like(xx.flatten())]).T
pts[:,2]+=0.05*np.sin(pts[:,0])+0.1*np.cos(pts[:,1]/2)# set the z value of the points
# second surface is a copy of first surface but with bigger z value 
pts2 = np.copy(pts)
pts2[:,2]+=1
pts2[:,3] = 1
pts = np.vstack([pts,pts2])
random_index = np.random.randint(0,len(pts),100)<jupyter_output><empty_output><jupyter_text>define fault surface as a straight line at x = 5<jupyter_code>fault = np.vstack([np.zeros(10)+5,np.linspace(0,10,10),np.zeros(10)-1,np.zeros(10)]).T<jupyter_output><empty_output><jupyter_text>Add these points into a pandas dataframe using the columns:* x* y* z* gx - gradient constraint x component* gy - gradient constraint y component* gz - gradient constraint z component* val - implicit function value* feature_name - identifying name for the geological feature<jupyter_code>from LoopStructural.utils.helper import empty_dataframe
import pandas as pd
data = pd.DataFrame(columns=['X','Y','Z','val','feature_name','gx','gy','gz'])#'Xempty_dataframe()
data = data.astype('float')
upper_df = pd.DataFrame(pts[random_index,:],columns=['X','Y','Z','val'])
upper_df['feature_name'] = 'upper'

fault_df = pd.DataFrame(fault,columns=['X','Y','Z','val'])
fault_df['gx'] = np.nan
fault_df['gy'] = np.nan
fault_df['gz'] = np.nan
fault_df.loc[len(fault_df),['X','Y','Z','gx','gy','gz']] = [5.,5.,-1.,0.7,0.,0.7]
fault_df['feature_name']='fault'

data = pd.concat([upper_df,fault_df])
data.reset_index(inplace=True)
data.loc[len(data),['X','Y','Z','gx','gy','gz','val','feature_name']] = [5.,5.,-2,.0,0,1,0,'lower']
origin = [0,0,-5]
maximum = [10,10,5]
model = GeologicalModel(origin,maximum,rescale=False)
model.set_model_data(data)
model.create_and_add_foliation('upper',interpolator_type='FDI',
                           nelements=1e4)
model.add_unconformity(model['upper'],value=0)
model.create_and_add_fault('fault',
                           1,
                           interpolatortype='FDI',
                           nelements=1e5,
                          fault_slip_vector=np.array([0.7,0.,-0.7]),
                          fault_center=np.array([5.,5.,-1.]),
                          fault_extent=1.,
                          fault_influence=1.,
                          fault_vertical_radius=1.,
                          buffer=0.5,
                           solver='pyamg'
                          )
model.create_and_add_foliation('lower',interpolator_type='FDI',
                           nelements=1e4)
model.update()<jupyter_output><empty_output><jupyter_text>Visualisation<jupyter_code>view = LavaVuModelViewer(model)
view.add_isosurface(model['upper'],value=0)
view.add_isosurface(model['lower'],slices=[0,-1,-2])

view.add_isosurface(model['upper'],value=1)
view.add_isosurface(model['fault'],value=0)

# view.add_vector_field(model['fault'][1],locations=model.regular_grid()[::100])
view.add_data(model['fault'][0])
view.add_data(model['fault'][1])
view.add_data(model['fault'][2])

view.add_data(model['upper'])

view.interactive()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10396409872849663, 'fraction_numerical': 0.0184848808633401, 'mean_word_length': 5.16600790513834, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '19850', 'license': 'CC0-1.0', 'n_tokens_mistral': 2947, 'n_tokens_neox': 2783, 'n_words': 884, 'path': 'notebooks/Building a geological model with LoopStructural.ipynb', 'repo_name': 'Loop3D/T21', 'size': '9359'}