docs: fixing new documentation so that data can be loaded

LoopStructural/datasets/__init__.py

@@ -17,3 +17,4 @@
 from ._base import load_unconformity
 from ._base import load_duplex
 from ._base import load_tabular_intrusion
+from ._base import load_geological_map_data

LoopStructural/datasets/_base.py

@@ -157,3 +157,69 @@ def load_tabular_intrusion():
     data = pd.read_csv(join(module_path, Path("data/tabular_intrusion.csv")))
     bb = np.array([[0, 0, 0], [5, 5, 5]])
     return data, bb
+
+
+def load_geological_map_data():
+    module_path = dirname(__file__)
+    contacts = pd.read_csv(
+        join(module_path, Path("data/geological_map_data/contacts.csv"))
+    )
+    stratigraphic_orientations = pd.read_csv(
+        join(
+            module_path, Path("data/geological_map_data/stratigraphic_orientations.csv")
+        )
+    )
+    stratigraphic_thickness = pd.read_csv(
+        join(module_path, Path("data/geological_map_data/stratigraphic_thickness.csv")),
+        skiprows=1,
+        names=["name", "thickness"],
+    )
+    stratigraphic_order = pd.read_csv(
+        join(module_path, Path("data/geological_map_data/stratigraphic_order.csv")),
+        skiprows=1,
+        names=["name", "order"],
+    )
+    bbox = pd.read_csv(
+        join(module_path, Path("data/geological_map_data/bbox.csv")),
+        index_col=0,
+        header=None,
+        names=["X", "Y", "Z"],
+    )
+    fault_properties = pd.read_csv(
+        join(module_path, Path("data/geological_map_data/fault_displacement.csv")),
+        index_col=0,
+    )
+    fault_edges = []
+    with open(
+        join(module_path, Path("data/geological_map_data/fault_edges.txt")), "r"
+    ) as f:
+        for l in f.read().split("\n"):
+            faults = l.split(",")
+            if len(faults) == 2:
+                fault_edges.append((faults[0], faults[1]))
+    fault_locations = pd.read_csv(
+        join(module_path, Path("data/geological_map_data/fault_locations.csv"))
+    )
+    fault_orientations = pd.read_csv(
+        join(module_path, Path("data/geological_map_data/fault_orientations.csv"))
+    )
+    return (
+        contacts,
+        stratigraphic_orientations,
+        stratigraphic_thickness,
+        stratigraphic_order,
+        bbox,
+        fault_locations,
+        fault_orientations,
+        fault_properties,
+        fault_edges,
+    )
+
+
+def load_fault_trace():
+    import geopandas
+
+    fault_trace = geopandas.read_file(
+        join(module_path, Path("data/fault_trace/fault_trace.shp"))
+    )
+    return fault_trace

examples/3_fault/fault_network.py

@@ -0,0 +1,125 @@
+"""
+3b. 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 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**| # | --------------- |-----| ------| -------|#  | ... | . | . | .
+
+try:
+    fault_trace = load_fault_trace()
+except:
+    print("Could not load fault trace")
+    exit()
+
+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() * 0.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 *= 0.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()