|
1 | 1 | """ |
2 | | -============================================ |
3 | | -Modelling a fault network in LoopStructural |
4 | | -
|
5 | | -============================================ |
| 2 | +3b. Modelling a fault network in LoopStructural |
| 3 | +=============================================== |
6 | 4 | Uses GeologicalModel, ProcessInputData and LavaVuModelViewer from LoopStructural library. Also using geopandas to read a shapefile, pandas, matplotlib and numpy.""" |
7 | 5 |
|
8 | 6 | import LoopStructural |
9 | 7 | LoopStructural.__version__ |
| 8 | + |
10 | 9 | from LoopStructural import GeologicalModel |
11 | 10 | from LoopStructural.modelling import ProcessInputData |
12 | 11 | from LoopStructural.visualisation import LavaVuModelViewer |
13 | 12 | import geopandas |
14 | 13 | import pandas as pd |
15 | 14 | import matplotlib.pyplot as plt |
16 | 15 | import numpy as np |
17 | | -# ~~~~~~~~~~~~~~~ |
| 16 | + |
| 17 | +############################## |
18 | 18 | # Read shapefile |
19 | | - |
20 | | -# ~~~~~~~~~~~~~~~ |
21 | | -# |
22 | | -# Read the shapefile and create a point for each node of the line |
23 | | -# |
24 | | -# | **fault_name** | **X** | **Y** | **Z**| |
25 | | -# | --------------- |-----| ------| -------| |
26 | | -# | ... | . | . | . |
| 19 | +# ~~~~~~~~~~~~~~ |
| 20 | +# # Read the shapefile and create a point for each node of the line # # | **fault_name** | **X** | **Y** | **Z**| # | --------------- |-----| ------| -------|# | ... | . | . | . |
| 21 | + |
27 | 22 | fault_trace = geopandas.read_file('fault_trace.shp') |
| 23 | + |
28 | 24 | faults = [] |
29 | 25 | for i in range(len(fault_trace)): |
30 | 26 | for x,y in zip(fault_trace.loc[i,:].geometry.xy[0],fault_trace.loc[i,:].geometry.xy[1]): |
31 | 27 | faults.append([fault_trace.loc[i,'fault_name'], x,y,np.random.random()*.4]) # better results if points aren't from a single plane |
32 | 28 | df = pd.DataFrame(faults,columns=['fault_name','X','Y','Z']) |
| 29 | + |
33 | 30 | fig, ax = plt.subplots() |
34 | 31 | ax.scatter(df['X'],df['Y']) |
35 | 32 | ax.axis('square') |
| 33 | + |
36 | 34 | scale = np.min([df['X'].max()-df['X'].min(),df['Y'].max()-df['Y'].min()]) |
37 | 35 | df['X']/=scale |
38 | 36 | df['Y']/=scale |
39 | 37 |
|
40 | | -# ~~~~~~~~~~~~~~~~~ |
| 38 | + |
| 39 | +############################## |
41 | 40 | # Orientation data |
42 | | - |
43 | | -# ~~~~~~~~~~~~~~~~~ |
| 41 | +# ~~~~~~~~~~~~~~~~ |
44 | 42 | # We can generate vertical dip data at the centre of the fault. |
| 43 | + |
45 | 44 | ori = [] |
46 | 45 | for f in df['fault_name'].unique(): |
47 | 46 | centre = df.loc[df['fault_name']==f,['X','Y','Z']].mean().to_numpy().tolist() |
|
51 | 50 | ori.append([f,*centre,*norm])#.extend(centre.extend(norm.tolist()))) |
52 | 51 | # fault_orientations = pd.DataFrame([[ |
53 | 52 | ori = pd.DataFrame(ori,columns=['fault_name','X','Y','Z','gx','gy','gz']) |
54 | | -# ~~~~~~~~~~~~~ |
| 53 | + |
| 54 | +############################## |
55 | 55 | # Model extent |
56 | | - |
57 | | -# ~~~~~~~~~~~~~ |
58 | | -# |
59 | | -# 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. |
| 56 | +# ~~~~~~~~~~~~ |
| 57 | +# # 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. |
| 58 | + |
60 | 59 | z = np.max([df['X'].max(),df['Y'].max()])- np.min([df['X'].min(),df['Y'].min()]) |
61 | 60 | z*=.2 |
62 | 61 | origin = [df['X'].min()-z,df['Y'].min()-z,-z] |
63 | 62 | maximum = [df['X'].max()+z,df['Y'].max()+z,z] |
64 | | -# ~~~~~~~~~~~~~~~~~~~~ |
| 63 | + |
| 64 | +############################## |
65 | 65 | # Setting up the data |
66 | | - |
67 | | -# ~~~~~~~~~~~~~~~~~~~~ |
68 | | -# The `ProcessInputData` class is used to convert common geological map components to the datastructures required by LoopStructural. |
69 | | -# |
70 | | -# To build a fault network we need to provide: |
71 | | -# * fault locations - a table of x,y,z, and the fault name |
72 | | -# * fault orientations - a table recording the orientation observations of the fault, e.g. strike, dip or normal vector and x,y,z, fault_name |
73 | | -# * origin - the origin of the model bounding box |
74 | | -# * maximum - the maximum extend of the model bounding box |
75 | | -# * fault_edges - list of intersection relationships between faults e.g. [('fault1','fault2')] indicates that there is a intersection between fault1 and fault2 |
76 | | -# * 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 |
77 | | -# * fault_properties (*optional*) - a pandas dataframe with any kwargs for the interpolator where the index is the fault name |
78 | | -# |
79 | | -# Below is an example of setting the number of interpolation elements for each fault |
80 | | -# ```Python |
81 | | -# fault_properties = pd.DataFrame([['fault_1',1e4], |
82 | | -# ['fault_2',1e4]], |
83 | | -# columns=['fault_name','nelements']).set_index('fault_name') |
84 | | -# ``` |
85 | | -# ~~~~~~~~~~~~~~~~~~~~~~~ |
| 66 | +# ~~~~~~~~~~~~~~~~~~~ |
| 67 | +# 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')# ``` |
| 68 | + |
| 69 | +############################## |
86 | 70 | # Modelling splay faults |
87 | | - |
88 | | -# ~~~~~~~~~~~~~~~~~~~~~~~ |
| 71 | +# ~~~~~~~~~~~~~~~~~~~~~~ |
89 | 72 | # 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$. |
| 73 | + |
90 | 74 | processor = ProcessInputData(fault_orientations=ori, |
91 | 75 | fault_locations=df, |
92 | 76 | origin=origin, |
93 | 77 | maximum=maximum, |
94 | 78 | fault_edges=[('fault_2','fault_1')], |
95 | 79 | fault_edge_properties=[{'angle':10}], |
96 | 80 | ) |
| 81 | + |
97 | 82 | model = GeologicalModel.from_processor(processor) |
98 | 83 | model.update() |
| 84 | + |
99 | 85 | view = LavaVuModelViewer(model) |
100 | 86 | for f in model.faults: |
101 | 87 | view.add_isosurface(f,slices=[0])# |
102 | 88 | view.rotation = [-50.92916488647461, -30.319700241088867, -20.521053314208984] |
103 | 89 | view.display() |
104 | | -# ~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 90 | + |
| 91 | +############################## |
105 | 92 | # Modelling abutting faults |
106 | | - |
107 | | -# ~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 93 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~ |
108 | 94 | # 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. |
| 95 | + |
109 | 96 | processor = ProcessInputData(fault_orientations=ori, |
110 | 97 | fault_locations=df, |
111 | 98 | origin=origin, |
112 | 99 | maximum=maximum, |
113 | 100 | fault_edges=[('fault_2','fault_1')], |
114 | 101 | fault_edge_properties=[{'angle':40}], |
115 | 102 | ) |
| 103 | + |
116 | 104 | model = GeologicalModel.from_processor(processor) |
| 105 | + |
117 | 106 | view = LavaVuModelViewer(model) |
118 | 107 | for f in model.faults: |
119 | 108 | view.add_isosurface(f,slices=[0])# |
|
0 commit comments