docs: adding fault network example

examples/3_fault/fault_trace.cpg

@@ -0,0 +1 @@
+UTF-8

examples/3_fault/fault_trace.prj

@@ -0,0 +1 @@
+PROJCS["GDA2020_MGA_Zone_53",GEOGCS["GCS_GDA2020",DATUM["GDA2020",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",135.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

examples/3_fault/plot_fault_network.py

@@ -0,0 +1,123 @@
+"""
+============================================
+Modelling a fault network in LoopStructural
+
+============================================
+Uses GeologicalModel, ProcessInputData and LavaVuModelViewer from LoopStructural library. Also using geopandas to read a shapefile, pandas, matplotlib and numpy."""
+
+import LoopStructural
+LoopStructural.__version__
+from LoopStructural import GeologicalModel
+from LoopStructural.modelling import ProcessInputData
+from LoopStructural.visualisation import LavaVuModelViewer
+import geopandas
+import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+# ~~~~~~~~~~~~~~~
+# Read shapefile
+
+# ~~~~~~~~~~~~~~~
+# 
+# Read the shapefile and create a point for each node of the line 
+# 
+# | **fault_name** | **X** | **Y** | **Z**| 
+# | --------------- |-----| ------| -------|
+#  | ... | . | . | .
+fault_trace = geopandas.read_file('fault_trace.shp')
+faults = []
+for i in range(len(fault_trace)):
+    for x,y in  zip(fault_trace.loc[i,:].geometry.xy[0],fault_trace.loc[i,:].geometry.xy[1]):
+        faults.append([fault_trace.loc[i,'fault_name'], x,y,np.random.random()*.4]) # better results if points aren't from a single plane
+df = pd.DataFrame(faults,columns=['fault_name','X','Y','Z'])
+fig, ax = plt.subplots()
+ax.scatter(df['X'],df['Y'])
+ax.axis('square')
+scale = np.min([df['X'].max()-df['X'].min(),df['Y'].max()-df['Y'].min()])
+df['X']/=scale
+df['Y']/=scale
+
+# ~~~~~~~~~~~~~~~~~
+# Orientation data
+
+# ~~~~~~~~~~~~~~~~~
+# We can generate vertical dip data at the centre of the fault.
+ori = []
+for f in df['fault_name'].unique():
+    centre = df.loc[df['fault_name']==f,['X','Y','Z']].mean().to_numpy().tolist()
+    tangent = df.loc[df['fault_name']==f,['X','Y','Z']].to_numpy()[0,:]-df.loc[df['fault_name']==f,['X','Y','Z']].to_numpy()[-1,:]
+    norm = tangent /np.linalg.norm(tangent) 
+    norm = norm.dot(np.array([[0,-1,0],[1,0,0],[0,0,0]]))
+    ori.append([f,*centre,*norm])#.extend(centre.extend(norm.tolist())))
+# fault_orientations = pd.DataFrame([[
+ori = pd.DataFrame(ori,columns=['fault_name','X','Y','Z','gx','gy','gz'])
+# ~~~~~~~~~~~~~
+# Model extent
+
+# ~~~~~~~~~~~~~
+# 
+# Calculate the bounding box for the model using the extent of the shapefiles. We make the Z coordinate 10% of the maximum x/y length.
+z = np.max([df['X'].max(),df['Y'].max()])- np.min([df['X'].min(),df['Y'].min()])
+z*=.2
+origin = [df['X'].min()-z,df['Y'].min()-z,-z]
+maximum = [df['X'].max()+z,df['Y'].max()+z,z]
+# ~~~~~~~~~~~~~~~~~~~~
+# Setting up the data
+
+# ~~~~~~~~~~~~~~~~~~~~
+# The `ProcessInputData` class is used to convert common geological map components to the datastructures required by LoopStructural.
+# 
+# To build a fault network we need to provide:
+# * fault locations - a table of x,y,z, and the fault name
+# * fault orientations - a table recording the orientation observations of the fault, e.g. strike, dip or normal vector and x,y,z, fault_name
+# * origin - the origin of the model bounding box
+# * maximum - the maximum extend of the model bounding box
+# * fault_edges - list of intersection relationships between faults e.g. [('fault1','fault2')] indicates that there is a intersection between fault1 and fault2
+# * fault_edge_properties - list of properties for the fault edges - this can be the type of intersection e.g. 'splay' or 'abut' or just the angle between the faults
+# * fault_properties (*optional*)  - a pandas dataframe with any kwargs for the interpolator where the index is the fault name 
+# 
+# Below is an example of setting the number of interpolation elements for each fault
+# ```Python
+# fault_properties = pd.DataFrame([['fault_1',1e4],
+#                                  ['fault_2',1e4]],
+#                                 columns=['fault_name','nelements']).set_index('fault_name')
+# ```
+# ~~~~~~~~~~~~~~~~~~~~~~~
+# Modelling splay faults
+
+# ~~~~~~~~~~~~~~~~~~~~~~~
+# A splay fault relationship is defined for any fault where the angle between the faults is less than $30^\circ$. In this example we specify the angle between the faults as $10^\circ$.
+processor = ProcessInputData(fault_orientations=ori, 
+                             fault_locations=df,
+                             origin=origin,
+                             maximum=maximum,
+                             fault_edges=[('fault_2','fault_1')],
+                             fault_edge_properties=[{'angle':10}],
+                            )
+model = GeologicalModel.from_processor(processor)
+model.update()
+view = LavaVuModelViewer(model)
+for f in model.faults:
+    view.add_isosurface(f,slices=[0])#
+view.rotation = [-50.92916488647461, -30.319700241088867, -20.521053314208984]
+view.display()
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Modelling abutting faults
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# In this exampe we will use the same faults but specify the angle between the faults as $40^\circ$ which will change the fault relationship to be abutting rather than splay. 
+processor = ProcessInputData(fault_orientations=ori, 
+                             fault_locations=df,
+                             origin=origin,
+                             maximum=maximum,
+                             fault_edges=[('fault_2','fault_1')],
+                             fault_edge_properties=[{'angle':40}],
+                            )
+model = GeologicalModel.from_processor(processor)
+view = LavaVuModelViewer(model)
+for f in model.faults:
+    view.add_isosurface(f,slices=[0])#
+    view.add_data(f[0],vectors=True)
+
+view.rotation = [-50.92916488647461, -30.319700241088867, -20.521053314208984]
+view.display()