1
1
"""
2
- Model 1 - Horizontal stratigraphy
3
- ==================================
2
+ Model 1 - Horizontal Stratigraphy
3
+ ===================================
4
4
5
- A simple example with horizontal layers
5
+ This example demonstrates how to build a basic geological model with horizontally stacked layers using GemPy,
6
+ an open-source Python library for implicit geological modeling. In this script, we will:
6
7
7
- This script demonstrates how to create a basic model with horizontally stacked layers using GemPy,
8
- a Python-based, open-source library for implicit geological modeling.
9
- """
8
+ - Set up and compute a simple geological model using input data (orientations and surface points).
9
+ - Visualize the model in 2D and 3D.
10
+ - Export model components (formations, orientations, surface points, and the regular grid) to VTK files.
11
+ - Extract mesh vertices and normals, save them to Excel files, and convert the vertex data to an XYZ file.
10
12
13
+ Each section includes detailed comments to explain its purpose and functionality.
14
+ """
11
15
12
16
# %%
13
17
# Import necessary libraries
14
18
import gempy as gp
15
19
import gempy_viewer as gpv
16
-
20
+ import pyvista as pv
21
+ import vtk
22
+ import pandas as pd
17
23
18
24
# sphinx_gallery_thumbnail_number = 2
19
25
20
26
# %%
21
- # Generate the model
22
- # Define the path to data
27
+ # Generate the Geological Model
28
+ # -----------------------------
29
+ # In this section, we define the data source and create a GeoModel instance.
30
+ # The model is defined with a specific spatial extent and resolution (refinement).
31
+ # Input data (orientations and surface points) are loaded from a remote repository.
32
+
33
+ # Define the path to the dataset hosted online
23
34
data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'
24
- # Create a GeoModel instance
35
+
36
+ # Create a GeoModel instance with the given project name, extent, and refinement level.
25
37
data = gp .create_geomodel (
26
38
project_name = 'horizontal' ,
27
39
extent = [0 , 1000 , 0 , 1000 , 0 , 1000 ],
31
43
path_to_surface_points = data_path + "/data/input_data/jan_models/model1_surface_points.csv"
32
44
)
33
45
)
34
- # Map geological series to surfaces
46
+
47
+ # Map a geological series to the corresponding surfaces.
48
+ # Here, the "Strat_Series" is associated with two formations ('rock2' and 'rock1'),
49
+ # which defines the stacking order of horizontal layers.
35
50
gp .map_stack_to_surfaces (
36
51
gempy_model = data ,
37
52
mapping_object = {"Strat_Series" : ('rock2' , 'rock1' )}
38
53
)
39
- # Compute the geological model
54
+
55
+ # Compute the geological model using the provided data and mappings.
40
56
gp .compute_model (data )
41
57
geo_data = data
42
58
43
59
# %%
44
- # Plot the initial geological model in the y direction without results
60
+ # 2D Visualization of the Geological Model
61
+ # -----------------------------------------
62
+ # The following plots provide 2D views of the model.
63
+ # - The first plot shows the initial model without computed results.
64
+ # - The subsequent plots display the computed geological data in the x and y directions.
65
+ # Note: The boundaries are hidden in the latter plots for a cleaner visualization.
66
+
67
+ # Plot the initial model in the y direction (without displaying computed results)
45
68
gpv .plot_2d (geo_data , direction = ['y' ], show_results = False )
46
69
47
- # Plot the result of the model in the x and y direction with data and without boundaries
70
+ # Plot the computed model results in the x and y directions, including data but excluding boundaries
48
71
gpv .plot_2d (geo_data , direction = ['x' ], show_data = True , show_boundaries = False )
49
72
gpv .plot_2d (geo_data , direction = ['y' ], show_data = True , show_boundaries = False )
50
73
51
74
# %%
52
- # Export to VTK model
53
- p = gpv .plot_3d (geo_data , show_data = True , show_results = True , show_boundaries = True ,image = True )
54
- p .surface_poly ['rock1' ].save ('rock1.vtk' ) # Save the vtk file for formation 1
55
- p .surface_poly ['rock2' ].save ('rock2.vtk' ) # Save the vtk file for formation 2
56
- p .orientations_mesh .save ('orientations.vtk' ) # Save the vtk file for the orientations
57
- p .surface_points_mesh .save ('surface_points.vtk' ) # Save the vtk file for the surface points
58
- box = p .regular_grid_actor .GetMapper ().GetInput () # Get the vtk file for the regular grid
75
+ # 3D Visualization and Export to VTK
76
+ # ----------------------------------
77
+ # Visualize the model in 3D with data, results, and boundaries.
78
+ # The generated 3D plot is also used to export various components of the model to VTK files.
79
+ p = gpv .plot_3d (geo_data , show_data = True , show_results = True , show_boundaries = True , image = True )
80
+
81
+ # Export different components of the model to VTK files:
82
+ p .surface_poly ['rock1' ].save ('rock1.vtk' ) # Save formation 'rock1'
83
+ p .surface_poly ['rock2' ].save ('rock2.vtk' ) # Save formation 'rock2'
84
+ p .orientations_mesh .save ('orientations.vtk' ) # Save the orientations mesh
85
+ p .surface_points_mesh .save ('surface_points.vtk' ) # Save the surface points mesh
86
+
87
+ # Retrieve and export the regular grid (volume representation) of the model
88
+ box = p .regular_grid_actor .GetMapper ().GetInput ()
59
89
box .save ('box.vtk' )
60
90
61
91
# %%
62
- import pyvista as pv
92
+ # Display the Exported VTK Files with PyVista
93
+ # -------------------------------------------
94
+ # Load and plot each exported VTK file to verify the export.
63
95
pv .read ('rock1.vtk' ).plot (show_edges = False )
64
96
pv .read ('rock2.vtk' ).plot (show_edges = False )
65
97
pv .read ('orientations.vtk' ).plot (show_edges = False )
66
98
pv .read ('surface_points.vtk' ).plot (show_edges = False )
67
99
pv .read ('box.vtk' ).plot (show_edges = False )
68
100
101
+
69
102
# %%
70
- # Export vertices of mesh
71
- from builtins import range
72
- import vtk
73
- import pandas as pd
103
+ # Extract and Save Mesh Vertices and Normals
104
+ # ------------------------------------------
105
+ # The following functions extract vertex coordinates and corresponding normal vectors
106
+ # from a mesh (vtkPolyData), and then save this data into Excel files for further analysis.
107
+
74
108
def generate_normals (polydata ):
109
+ """
110
+ Generate normals for the given polydata if they are not already computed.
111
+
112
+ Parameters:
113
+ polydata (vtk.vtkPolyData): Input polydata for which normals are to be computed.
114
+
115
+ Returns:
116
+ vtk.vtkPolyData: The polydata with computed point normals.
117
+ """
75
118
normal_generator = vtk .vtkPolyDataNormals ()
76
119
normal_generator .SetInputData (polydata )
77
120
normal_generator .ComputePointNormalsOn ()
78
121
normal_generator .ComputeCellNormalsOff ()
79
122
normal_generator .Update ()
80
123
return normal_generator .GetOutput ()
81
124
125
+
82
126
def get_vertices_and_normals (mesh ):
127
+ """
128
+ Extract vertices and normals from a mesh.
129
+
130
+ Parameters:
131
+ mesh: The mesh object (typically from a formation) to extract surface vertices and normals.
83
132
133
+ Returns:
134
+ tuple: Two lists containing vertices and normals respectively.
135
+ """
136
+ # Extract the surface of the mesh
84
137
surface_mesh = mesh .extract_surface ()
85
138
polydata = surface_mesh
86
139
87
- # Generate normals if not present
140
+ # Ensure normals are computed for the polydata
88
141
polydata_with_normals = generate_normals (polydata )
89
142
90
- # Get points (vertices)
143
+ # Extract points (vertices)
91
144
points = polydata_with_normals .GetPoints ()
92
- vertices = []
93
- for i in range (points .GetNumberOfPoints ()):
94
- vertices .append (points .GetPoint (i ))
145
+ vertices = [points .GetPoint (i ) for i in range (points .GetNumberOfPoints ())]
95
146
96
- # Get normals
147
+ # Extract normals associated with the points
97
148
normals_array = polydata_with_normals .GetPointData ().GetNormals ()
98
- normals = []
99
- for i in range (normals_array .GetNumberOfTuples ()):
100
- normals .append (normals_array .GetTuple (i ))
149
+ normals = [normals_array .GetTuple (i ) for i in range (normals_array .GetNumberOfTuples ())]
101
150
102
151
return vertices , normals
103
152
153
+
104
154
def save_to_excel (vertices , normals , vertices_file , normals_file ):
105
- # Create DataFrames
155
+ """
156
+ Save the vertices and normals to separate Excel files.
157
+
158
+ Parameters:
159
+ vertices (list): List of vertex coordinates.
160
+ normals (list): List of normal vectors.
161
+ vertices_file (str): File name for saving vertices.
162
+ normals_file (str): File name for saving normals.
163
+ """
164
+ # Create DataFrames from the lists
106
165
vertices_df = pd .DataFrame (vertices , columns = ['X' , 'Y' , 'Z' ])
107
166
normals_df = pd .DataFrame (normals , columns = ['x' , 'y' , 'z' ])
108
167
109
- # Save to Excel files
168
+ # Write the DataFrames to Excel files
110
169
vertices_df .to_excel (vertices_file , index = False )
111
- normals_df .to_excel
170
+ normals_df .to_excel (normals_file , index = False )
171
+
172
+
173
+ # Extract vertices and normals from the mesh of formation 'rock1'
112
174
mesh = p .surface_poly ['rock1' ]
113
175
vertices , normals = get_vertices_and_normals (mesh )
176
+
177
+ # Create DataFrames for later processing or verification
114
178
vertices_df = pd .DataFrame (vertices , columns = ['X' , 'Y' , 'Z' ])
115
179
normals_df = pd .DataFrame (normals , columns = ['x' , 'y' , 'z' ])
116
- # Save to Excel filesthe
180
+
181
+ # Define file names for the Excel outputs
117
182
vertices_file = "rock1_vertices.xlsx"
118
183
normals_file = "rock1_norms.xlsx"
184
+
185
+ # Save the extracted data to Excel files
119
186
save_to_excel (vertices , normals , vertices_file , normals_file )
187
+
188
+ # Optionally, read back the Excel files to verify their content
120
189
pd .read_excel (vertices_file )
121
190
pd .read_excel (normals_file )
122
191
192
+
123
193
# %%
124
- # Convert the DataFrame to an XYZ file
194
+ # Convert Vertices DataFrame to an XYZ File
195
+ # -----------------------------------------
196
+ # This function converts a DataFrame containing vertex coordinates to an XYZ format file,
197
+ # which is a simple text file with one point per line.
198
+
125
199
def dataframe_to_xyz (df , filename ):
200
+ """
201
+ Write vertex coordinates from a DataFrame to an XYZ file.
202
+
203
+ Parameters:
204
+ df (pandas.DataFrame): DataFrame containing 'X', 'Y', and 'Z' columns.
205
+ filename (str): Output filename for the XYZ file.
206
+ """
126
207
with open (filename , 'w' ) as f :
127
208
for index , row in df .iterrows ():
128
209
f .write (f"{ row ['X' ]} { row ['Y' ]} { row ['Z' ]} \n " )
129
210
130
- # Specify the filename
211
+
212
+ # Specify the output file name for the XYZ file
131
213
filename = "output.xyz"
132
214
133
- # Call the function to write the DataFrame to an XYZ file
134
- dataframe_to_xyz (vertices_df , filename )
215
+ # Convert the vertices DataFrame to an XYZ file
216
+ dataframe_to_xyz (vertices_df , filename )
0 commit comments