/
plot_2_surface_modelling.py
184 lines (152 loc) · 6.13 KB
/
plot_2_surface_modelling.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
"""
1b. Implicit surface modelling
===============================
This tutorial will demonstrate how to create an implicit surface
representation of surfaces from a combination of orientation and
location observations.
Implicit surface representation involves finding an unknown function
where :math:`f(x,y,z)` matches observations of the surface geometry. We
generate a scalar field where the scalar value is the distance away from
a reference horizon. The reference horizon is arbritary and can either
be:
- a single geological surface where the scalar field would represent
the signed distance away from this surface. (above the surface
positive and below negative)
- Where multiple conformable horizons are observed the same scalar
field can be used to represent these surfaces and the thickness of
the layers is used to determine the relative scalar value for each
surface
This tutorial will demonstrate both of these approaches for modelling a
number of horizons picked from seismic data sets, by following the next
steps: 1. Creation of a geological model, which includes: \*
Presentation and visualization of the data \* Addition of a geological
feature, which in this case is the stratigraphy of the model. 2.
Visualization of the scalar field.
"""
#########################################################################
# Imports
# ~~~~~~~
# Import the required objects from LoopStructural for visualisation and
# model building
from LoopStructural import GeologicalModel
from LoopStructural.visualisation import LavaVuModelViewer
from LoopStructural.datasets import load_claudius # demo data
import numpy as np
######################################################################
# The data for this example can be imported from the example datasets
# module in loopstructural.
#
data, bb = load_claudius()
data["val"].unique()
######################################################################
# GeologicalModel
# ~~~~~~~~~~~~~~~
#
# To create a ``GeologicalModel`` the origin (lower left) and max extent
# (upper right) of the model area need to be specified. In this example
# these are provided in the bb array.
#
model = GeologicalModel(bb[0, :], bb[1, :])
######################################################################
# A pandas dataframe with appropriate columns can be used to link the data
# to the geological model.
#
# * ``X`` is the x coordinate
# * ``Y`` is the y # coordinate
# * ``Z`` is the z coordinate
# * ``feature_name`` is a name to link the data to a model object
# * ``val`` is the value of the scalar field which represents the
# distance from a reference horizon. It is comparable
# to the relative thickness
#
# * ``nx`` is the x component of the normal vector to the surface gradient
# * ``ny`` is the y component of the normal vector to the surface gradient
# * ``nz`` is the z component of the normal vector to the surface gradeint
# * ``strike`` is the strike angle
# * ``dip`` is the dip angle
#
# Having a look at the data for this example by looking at the top of the
# dataframe and then using a 3D plot
#
data["feature_name"].unique()
viewer = LavaVuModelViewer(background="white")
viewer.add_value_data(
data[~np.isnan(data["val"])][["X", "Y", "Z"]],
data[~np.isnan(data["val"])]["val"],
name="value points",
)
viewer.add_vector_data(
data[~np.isnan(data["nx"])][["X", "Y", "Z"]],
data[~np.isnan(data["nx"])][["nx", "ny", "nz"]],
name="orientation points",
)
viewer.rotate([-85.18760681152344, 42.93233871459961, 0.8641873002052307])
viewer.display()
######################################################################
# The pandas dataframe can be linked to the ``GeologicalModel`` using
# ``.set_model_data(dataframe``
#
model.set_model_data(data)
######################################################################
# The ``GeologicalModel`` contains an ordered list of the different
# geological features within a model and how these features interact. This
# controls the topology of the different geological features in the model.
# Different geological features can be added to the geological model such
# as:
# * Foliations
# * Faults
# * Unconformities
# * Folded foliations
# * Structural Frames
#
# In this example we will only add a foliation using the function
#
# .. code:: python
#
# model.create_and_add_foliation(name)
#
# where name is the name in the ``feature_name`` field, other parameters we
# specified are the:
# * ``interpolatortype`` - we can either use a
# PiecewiseLinearInterpolator ``PLI``, a FiniteDifferenceInterpolator
# ``FDI`` or a radial basis interpolator ``surfe``
# * ``nelements - int`` is the how many elements are used to discretize the resulting solution
# * ``buffer - float`` buffer percentage around the model area
# * ``solver`` - the algorithm to solve the least squares problem e.g.
# ``lu`` for lower upper decomposition, ``cg`` for conjugate gradient,
# ``pyamg`` for an algorithmic multigrid solver
# * ``damp - bool`` - whether to add a small number to the diagonal of the interpolation
# matrix for discrete interpolators - this can help speed up the solver
# and makes the solution more stable for some interpolators
#
vals = [0, 60, 250, 330, 600]
strat_column = {"strati": {}}
for i in range(len(vals) - 1):
strat_column["strati"]["unit_{}".format(i)] = {
"min": vals[i],
"max": vals[i + 1],
"id": i,
}
model.set_stratigraphic_column(strat_column)
strati = model.create_and_add_foliation(
"strati",
interpolatortype="FDI", # try changing this to 'PLI'
nelements=1e4, # try changing between 1e3 and 5e4
buffer=0.3,
solver="pyamg",
damp=True,
)
######################################################################
# Plot the surfaces
# ------------------------------------
viewer = LavaVuModelViewer(model)
viewer.add_model_surfaces(cmap="tab20")
viewer.rotate([-85.18760681152344, 42.93233871459961, 0.8641873002052307])
viewer.display()
######################################################################
# Plot block diagram
# -------------------
viewer = LavaVuModelViewer(model)
viewer.add_model(cmap="tab20")
viewer.rotate([-85.18760681152344, 42.93233871459961, 0.8641873002052307])
viewer.display()