Skip to content

Commit bcc07ad

Browse files
LeguarkLeguark
committed
[ENH] Adding function to create elements from subsurface wells easier
1 parent a051c25 commit bcc07ad

File tree

3 files changed

+72
-75
lines changed

3 files changed

+72
-75
lines changed

‎gempy/API/__init__.py

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
11
# Initialization API
22
from .initialization_API import (
33
create_data_legacy,
4-
create_geomodel
4+
create_geomodel,
5+
structural_elements_from_borehole_set
56
)
67

78
# Compute API
@@ -62,12 +63,13 @@
6263
from gempy_engine.modules.geophysics.gravity_gradient import calculate_gravity_gradient
6364

6465
__all__ = [
65-
'create_data_legacy', 'create_geomodel', 'compute_model', 'compute_model_at', 'map_stack_to_surfaces',
66-
'set_section_grid', 'set_active_grid', 'set_topography_from_random', 'set_topography_from_file', 'set_topography_from_subsurface_structured_grid', 'set_topography_from_arrays',
67-
'set_custom_grid', 'set_centered_grid',
68-
'generate_example_model', 'set_fault_relation', 'set_is_fault', 'set_is_finite_fault',
69-
'add_surface_points', 'add_orientations', 'delete_surface_points', 'delete_orientations',
70-
'create_orientations_from_surface_points_coords', 'modify_surface_points', 'modify_orientations',
71-
'add_structural_group', 'remove_structural_group_by_index', 'remove_structural_group_by_name', 'remove_element_by_name',
72-
'calculate_gravity_gradient'
66+
'create_data_legacy', 'create_geomodel', 'structural_elements_from_borehole_set',
67+
'compute_model', 'compute_model_at', 'map_stack_to_surfaces',
68+
'set_section_grid', 'set_active_grid', 'set_topography_from_random', 'set_topography_from_file', 'set_topography_from_subsurface_structured_grid', 'set_topography_from_arrays',
69+
'set_custom_grid', 'set_centered_grid',
70+
'generate_example_model', 'set_fault_relation', 'set_is_fault', 'set_is_finite_fault',
71+
'add_surface_points', 'add_orientations', 'delete_surface_points', 'delete_orientations',
72+
'create_orientations_from_surface_points_coords', 'modify_surface_points', 'modify_orientations',
73+
'add_structural_group', 'remove_structural_group_by_index', 'remove_structural_group_by_name', 'remove_element_by_name',
74+
'calculate_gravity_gradient'
7375
]

‎gempy/API/initialization_API.py

Lines changed: 45 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,13 @@
11
import warnings
2-
from typing import Union
2+
from typing import Union, Hashable
33

44
import numpy as np
55
from numpy import ndarray
66

77
from gempy.API.io_API import read_surface_points, read_orientations
88
from gempy_engine.core.data import InterpolationOptions
9+
from ..optional_dependencies import require_subsurface
10+
from ..core.data import StructuralElement
911
from ..core.data.geo_model import GeoModel
1012
from ..core.data.grid import Grid
1113
from ..core.data.importer_helper import ImporterHelper
@@ -78,6 +80,48 @@ def create_geomodel(
7880
return geo_model
7981

8082

83+
def structural_elements_from_borehole_set(
84+
borehole_set: "subsurface.core.geological_formats.BoreholeSet",
85+
elements_dict: dict) -> list[StructuralElement]:
86+
"""
87+
Creates a list of StructuralElement instances from a borehole set and a dictionary of elements.
88+
89+
Args:
90+
borehole_set (subsurface.core.geological_formats.boreholes.boreholes): The borehole set.
91+
elements_dict (dict): A dictionary of elements.
92+
93+
Returns:
94+
list[StructuralElement]: A list of StructuralElement instances.
95+
"""
96+
97+
ss = require_subsurface()
98+
borehole_set: ss.core.geological_formats.BoreholeSet
99+
100+
elements = []
101+
component_lith: dict[Hashable, np.ndarray] = borehole_set.get_top_coords_for_each_lith()
102+
103+
for name, properties in elements_dict.items():
104+
top_coordinates = component_lith.get(properties['top_lith'])
105+
if top_coordinates is None:
106+
raise ValueError(f"Top lithology {properties['top_lith']} not found in borehole set.")
107+
108+
element = StructuralElement(
109+
name=name,
110+
id=properties['id'],
111+
color=properties['color'],
112+
surface_points=SurfacePointsTable.from_arrays(
113+
x=top_coordinates[:, 0],
114+
y=top_coordinates[:, 1],
115+
z=top_coordinates[:, 2],
116+
names=[name],
117+
name_id_map={name: properties['id']}
118+
),
119+
orientations=OrientationsTable(np.zeros(0, dtype=OrientationsTable.dt))
120+
)
121+
elements.append(element)
122+
return elements
123+
124+
81125
def create_data_legacy(
82126
*,
83127
project_name: str = 'default_project',

‎test/test_modules/test_pile/test_stratigraphic_pile.py

Lines changed: 16 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -75,74 +75,22 @@ def test_structural_elements(self, borehole_set: BoreholeSet):
7575
)
7676
pv_plot([s], image_2d=False, cmap="tab20c")
7777

78-
vertex_attributes: pd.DataFrame = borehole_trajectory.data.points_attributes
79-
unique_lith_codes = vertex_attributes['component lith'].unique()
80-
81-
component_lith = borehole_set.compute_tops()
82-
83-
pleistozen = gp.data.StructuralElement(
84-
name="Pleistozen",
85-
id=10_000,
86-
color="#f9f97f",
87-
surface_points=gp.data.SurfacePointsTable(np.empty(0, dtype=gp.data.SurfacePointsTable.dt)),
88-
orientations=gp.data.OrientationsTable(np.zeros(0, dtype=gp.data.OrientationsTable.dt))
89-
)
90-
91-
kreide = gp.data.StructuralElement(
92-
name="Kreide",
93-
id=30_000,
94-
color="#a6d84a",
95-
surface_points=gp.data.SurfacePointsTable(np.empty(0, dtype=gp.data.SurfacePointsTable.dt)),
96-
orientations=gp.data.OrientationsTable(np.zeros(0, dtype=gp.data.OrientationsTable.dt))
97-
)
98-
99-
trias = gp.data.StructuralElement(
100-
name="Trias",
101-
id=50_000,
102-
color="#a4469f",
103-
surface_points=gp.data.SurfacePointsTable(np.empty(0, dtype=gp.data.SurfacePointsTable.dt)),
104-
orientations=gp.data.OrientationsTable(np.zeros(0, dtype=gp.data.OrientationsTable.dt))
105-
)
106-
107-
perm = gp.data.StructuralElement(
108-
name="Perm",
109-
id=60_000,
110-
color="#f4a142",
111-
surface_points=gp.data.SurfacePointsTable(np.empty(0, dtype=gp.data.SurfacePointsTable.dt)),
112-
orientations=gp.data.OrientationsTable(np.zeros(0, dtype=gp.data.OrientationsTable.dt))
113-
)
114-
115-
rotliegend_id = 62_000
116-
rotliegend_xyz = component_lith[rotliegend_id]
117-
118-
# Add the id
119-
rotliegend_surface_points = gp.data.SurfacePointsTable.from_arrays(
120-
x=rotliegend_xyz[:, 0],
121-
y=rotliegend_xyz[:, 1],
122-
z=rotliegend_xyz[:, 2],
123-
names=["Rotliegend"],
124-
name_id_map={"Rotliegend": rotliegend_id}
125-
)
126-
127-
rotliegend = gp.data.StructuralElement(
128-
name="Rotliegend",
129-
id=rotliegend_id,
130-
color="#bb825b",
131-
surface_points=rotliegend_surface_points,
132-
orientations=gp.data.OrientationsTable(np.zeros(0, dtype=gp.data.OrientationsTable.dt))
133-
)
134-
135-
devon = gp.data.StructuralElement(
136-
name="Devon",
137-
id=80_000,
138-
color="#969594",
139-
surface_points=gp.data.SurfacePointsTable(np.empty(0, dtype=gp.data.SurfacePointsTable.dt)),
140-
orientations=gp.data.OrientationsTable(np.zeros(0, dtype=gp.data.OrientationsTable.dt))
78+
elements = gp.structural_elements_from_borehole_set(
79+
borehole_set=borehole_set,
80+
elements_dict={
81+
# "Pleistozen": {"id": 10_000, "color": "#f9f97f", "top_lith": 10_000},
82+
# "Kreide": {"id": 30_000, "color": "#a6d84a", "top_lith": 30_000},
83+
# "Trias": {"id": 50_000, "color": "#a4469f", "top_lith": 50_000},
84+
# "Perm": {"id": 60_000, "color": "#f4a142", "top_lith": 60_000},
85+
"Rotliegend": {"id": 62_000, "color": "#bb825b", "top_lith": 62_000},
86+
# "Devon": {"id": 80_000, "color": "#969594", "top_lith": 80_000}
87+
}
14188
)
142-
89+
90+
14391
group = gp.data.StructuralGroup(
14492
name="Stratigraphic Pile",
145-
elements=[rotliegend],
93+
elements=elements,
14694
structural_relation=gp.data.StackRelationType.ERODE
14795
)
14896
structural_frame = gp.data.StructuralFrame(
@@ -151,6 +99,9 @@ def test_structural_elements(self, borehole_set: BoreholeSet):
15199
)
152100
print(group)
153101

102+
103+
component_lith = borehole_set.get_top_coords_for_each_lith()
104+
rotliegend_xyz = component_lith[62_000]
154105
extent_from_data = rotliegend_xyz.min(axis=0), rotliegend_xyz.max(axis=0)
155106

156107
geo_model = gp.data.GeoModel(

0 commit comments

Comments
 (0)