Skip to content

Commit eb4049b

Browse files
LeguarkLeguark
committed
[DOC] Adding all the code for vector-geology examples
1 parent 598351e commit eb4049b

File tree

7 files changed

+777
-0
lines changed

7 files changed

+777
-0
lines changed

‎.env.example

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
LIQUID_EARTH_API_TOKEN = "foobar-1499cbe49808c8cfef14303d2aaacXDx-w1vc1bD065Wj0ABi3HD2fzQ"
2+
BOREHOLES_SOUTH_FOLDER = "/mnt/d/OneDrive - .../boreholes_south"

‎examples/examples/README.rst

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
Examples
2+
========
3+
4+
This examples intend to be a gallery of what other people has done with GemPy.
5+
Hopefully with the help of the community we can create a complete library of
6+
structures that can serve as template for new users.

‎examples/examples/__init__.py

Whitespace-only changes.

‎examples/examples/_aux_func_gravity_inversion.py

Lines changed: 194 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,194 @@
1+
import numpy as np
2+
from sklearn.linear_model import LinearRegression
3+
4+
def calculate_scale_shift(a: np.ndarray, b: np.ndarray) -> tuple:
5+
# Reshape arrays for sklearn
6+
a_reshaped = a.reshape(-1, 1)
7+
b_reshaped = b.reshape(-1, 1)
8+
9+
# Linear regression
10+
model = LinearRegression().fit(a_reshaped, b_reshaped)
11+
12+
# Get scale and shift
13+
s = model.coef_[0][0]
14+
c = model.intercept_[0]
15+
16+
return s, c
17+
18+
19+
def gaussian_kernel(locations, length_scale, variance):
20+
import torch
21+
# Compute the squared Euclidean distance between each pair of points
22+
locations = torch.tensor(locations.values)
23+
distance_squared = torch.cdist(locations, locations, p=2).pow(2)
24+
# Compute the covariance matrix using the Gaussian kernel
25+
covariance_matrix = variance * torch.exp(-0.5 * distance_squared / length_scale**2)
26+
return covariance_matrix
27+
28+
29+
import numpy as np
30+
31+
import gempy as gp
32+
import xarray as xr
33+
from vector_geology.model_building_functions import optimize_nuggets_for_group
34+
35+
36+
def initialize_geo_model(structural_elements: list[gp.data.StructuralElement], extent: list[float],
37+
topography: xr.Dataset, load_nuggets: bool = False
38+
) -> gp.data.GeoModel:
39+
structural_group_red = gp.data.StructuralGroup(
40+
name="Red",
41+
elements=[structural_elements[i] for i in [0, 4, 8]],
42+
structural_relation=gp.data.StackRelationType.ERODE
43+
)
44+
45+
# Any, Probably we can decimize this an extra notch
46+
structural_group_green = gp.data.StructuralGroup(
47+
name="Green",
48+
elements=[structural_elements[i] for i in [5]],
49+
structural_relation=gp.data.StackRelationType.ERODE
50+
)
51+
52+
# Blue range 2 cov 4
53+
structural_group_blue = gp.data.StructuralGroup(
54+
name="Blue",
55+
elements=[structural_elements[i] for i in [2, 3]],
56+
structural_relation=gp.data.StackRelationType.ERODE
57+
)
58+
59+
structural_group_intrusion = gp.data.StructuralGroup(
60+
name="Intrusion",
61+
elements=[structural_elements[i] for i in [1]],
62+
structural_relation=gp.data.StackRelationType.ERODE
63+
)
64+
65+
structural_groups = [structural_group_intrusion, structural_group_green, structural_group_blue, structural_group_red]
66+
structural_frame = gp.data.StructuralFrame(
67+
structural_groups=structural_groups[2:],
68+
color_gen=gp.data.ColorsGenerator()
69+
)
70+
# TODO: If elements do not have color maybe loop them on structural frame constructor?
71+
72+
geo_model: gp.data.GeoModel = gp.create_geomodel(
73+
project_name='Tutorial_ch1_1_Basics',
74+
extent=extent,
75+
resolution=[20, 10, 20],
76+
refinement=5, # * Here we define the number of octree levels. If octree levels are defined, the resolution is ignored.
77+
structural_frame=structural_frame
78+
)
79+
80+
if topography is not None:
81+
gp.set_topography_from_arrays(
82+
grid=geo_model.grid,
83+
xyz_vertices=topography.vertex.values
84+
)
85+
86+
if load_nuggets:
87+
import os
88+
89+
project_root = os.getcwd()
90+
path_to_temp = os.path.join(project_root, "../temp")
91+
apply_optimized_nuggets(
92+
geo_model=geo_model,
93+
loaded_nuggets_red=(np.load(path_to_temp + "/nuggets_Red.npy")),
94+
loaded_nuggets_blue=(np.load(path_to_temp + "/nuggets_Blue.npy")),
95+
loaded_nuggets_green=(np.load(path_to_temp + "/nuggets_Green.npy"))
96+
)
97+
98+
geo_model.structural_frame.get_element_by_name("KKR").color = "#A46283"
99+
geo_model.structural_frame.get_element_by_name("LGR").color = "#6394A4"
100+
geo_model.structural_frame.get_element_by_name("WAL").color = "#72A473"
101+
geo_model.structural_frame.get_element_by_name("ABL").color = "#1D3943"
102+
geo_model.structural_frame.basement_color = "#8B4220"
103+
104+
geo_model.update_transform()
105+
106+
return geo_model
107+
108+
109+
def optimize_nuggets_for_whole_project(geo_model: gp.data.GeoModel):
110+
geo_model.interpolation_options.kernel_options.range = 0.7
111+
geo_model.interpolation_options.kernel_options.c_o = 4
112+
optimize_nuggets_for_group(
113+
geo_model=geo_model,
114+
structural_group=geo_model.structural_frame.get_group_by_name('Red'),
115+
plot_evaluation=False,
116+
plot_result=True
117+
)
118+
geo_model.interpolation_options.kernel_options.range = 2
119+
geo_model.interpolation_options.kernel_options.c_o = 4
120+
optimize_nuggets_for_group(
121+
geo_model=geo_model,
122+
structural_group=geo_model.structural_frame.get_group_by_name('Blue'),
123+
plot_evaluation=False,
124+
plot_result=False
125+
)
126+
if False:
127+
optimize_nuggets_for_group(
128+
geo_model=geo_model,
129+
structural_group=geo_model.structural_frame.get_group_by_name('Green'),
130+
plot_evaluation=False,
131+
plot_result=True
132+
)
133+
134+
135+
def apply_optimized_nuggets(geo_model: gp.data.GeoModel, loaded_nuggets_red, loaded_nuggets_blue, loaded_nuggets_green):
136+
gp.modify_surface_points(
137+
geo_model,
138+
slice=None,
139+
elements_names=[element.name for element in geo_model.structural_frame.get_group_by_name('Red').elements],
140+
nugget=loaded_nuggets_red
141+
)
142+
if True: # Ignore OB
143+
gp.modify_surface_points(
144+
geo_model,
145+
slice=None,
146+
elements_names=[element.name for element in geo_model.structural_frame.get_group_by_name('Blue').elements],
147+
nugget=loaded_nuggets_blue
148+
)
149+
if False:
150+
gp.modify_surface_points(
151+
geo_model,
152+
slice=None,
153+
elements_names=[element.name for element in geo_model.structural_frame.get_group_by_name('Green').elements],
154+
nugget=loaded_nuggets_green
155+
)
156+
157+
158+
def setup_geophysics(env_path: str, geo_model: gp.data.GeoModel):
159+
import pandas as pd
160+
from dotenv import dotenv_values
161+
config = dotenv_values()
162+
163+
df = pd.read_csv(
164+
filepath_or_buffer=config.get(env_path),
165+
sep=',',
166+
header=0
167+
)
168+
169+
# Remove the items that have X > 5650000
170+
df = df[df['X'] < 565000]
171+
172+
interesting_columns = df[['X', 'Y', 'Bouguer_267_complete']]
173+
# %%
174+
device_location = interesting_columns[['X', 'Y']]
175+
# stack 0 to the z axis
176+
device_location['Z'] = 0
177+
178+
gp.set_centered_grid(
179+
grid=geo_model.grid,
180+
centers=device_location,
181+
resolution=np.array([10, 10, 15]),
182+
radius=np.array([5000, 5000, 5000])
183+
)
184+
185+
gravity_gradient = gp.calculate_gravity_gradient(geo_model.grid.centered_grid)
186+
187+
# %%
188+
from gempy_engine.core.backend_tensor import BackendTensor
189+
geo_model.geophysics_input = gp.data.GeophysicsInput(
190+
tz=BackendTensor.t.array(gravity_gradient),
191+
densities=BackendTensor.t.array([2.61, 2.92, 3.1, 2.92, 2.61, 2.61]),
192+
)
193+
194+
return interesting_columns

‎examples/examples/_aux_func_petrophysics_inversion.py

Lines changed: 177 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,177 @@
1+
import numpy as np
2+
import os
3+
import pandas as pd
4+
import pyvista
5+
6+
import gempy as gp
7+
import gempy_viewer as gpv
8+
import subsurface as ss
9+
from subsurface.modules.visualization import init_plotter, to_pyvista_points, to_pyvista_line
10+
11+
12+
def generate_spremberg_model(
13+
borehole_set: ss.core.geological_formats.BoreholeSet,
14+
elements_to_gempy: dict[str, dict[str, str]], plot: bool = False) -> gp.data.GeoModel:
15+
elements: list[gp.data.StructuralElement] = gp.structural_elements_from_borehole_set(
16+
borehole_set=borehole_set,
17+
elements_dict=elements_to_gempy
18+
)
19+
20+
group = gp.data.StructuralGroup(
21+
name="Stratigraphic Pile",
22+
elements=elements,
23+
structural_relation=gp.data.StackRelationType.ERODE
24+
)
25+
structural_frame = gp.data.StructuralFrame(
26+
structural_groups=[group],
27+
color_gen=gp.data.ColorsGenerator()
28+
)
29+
30+
# %%
31+
# Get the extent from the borehole set
32+
all_surface_points_coords: gp.data.SurfacePointsTable = structural_frame.surface_points_copy
33+
extent_from_data = all_surface_points_coords.xyz.min(axis=0), all_surface_points_coords.xyz.max(axis=0)
34+
35+
# Calculate point_y_axis
36+
n_octree_levels = 5
37+
dz = extent_from_data[1][2] - extent_from_data[0][2]
38+
zmin = extent_from_data[0][2] - dz * 0.1
39+
zmax = extent_from_data[1][2] + dz * 0.1
40+
regular_grid = gp.data.grid.RegularGrid.from_corners_box(
41+
pivot=(5_478_256.5, 5_698_528.946534388),
42+
point_x_axis=((5_483_077.527386775, 5_710_030.2446156405)),
43+
distance_point3=35_000,
44+
zmin=zmin,
45+
zmax=zmax,
46+
resolution=np.array([2 ** n_octree_levels] * 3),
47+
plot=True
48+
)
49+
50+
interpolation_options = gp.data.InterpolationOptions(
51+
range=5,
52+
c_o=10,
53+
mesh_extraction=True,
54+
number_octree_levels=n_octree_levels,
55+
)
56+
57+
grid = gp.data.grid.Grid()
58+
grid.set_octree_grid(regular_grid, interpolation_options.evaluation_options)
59+
60+
geo_model = gp.data.GeoModel(
61+
name="Stratigraphic Pile",
62+
structural_frame=structural_frame,
63+
grid=grid,
64+
interpolation_options=interpolation_options,
65+
)
66+
67+
gempy_plot = gpv.plot_3d(
68+
model=geo_model,
69+
kwargs_pyvista_bounds={
70+
'show_xlabels': False,
71+
'show_ylabels': False,
72+
},
73+
show=True,
74+
image=True
75+
)
76+
77+
if plot:
78+
add_wells_plot(borehole_set, gempy_plot, geo_model.grid)
79+
80+
return geo_model
81+
82+
83+
def get_spremberg_borehole_set() -> ss.core.geological_formats.BoreholeSet:
84+
borehole_set: ss.core.geological_formats.BoreholeSet = _read_spremberg_borehole_set()
85+
return borehole_set
86+
87+
88+
def add_wells_plot(borehole_set, gempy_plot, grid: gp.data.grid.Grid):
89+
sp_mesh: pyvista.PolyData = gempy_plot.surface_points_mesh
90+
well_mesh = to_pyvista_line(
91+
line_set=borehole_set.combined_trajectory,
92+
active_scalar="lith_ids",
93+
radius=10
94+
)
95+
units_limit = [0, 20]
96+
collars = to_pyvista_points(
97+
borehole_set.collars.collar_loc,
98+
)
99+
pyvista_plotter = init_plotter(ve=10)
100+
pyvista_plotter.show_bounds(all_edges=False)
101+
pyvista_plotter.add_mesh(
102+
well_mesh.threshold(units_limit),
103+
cmap="tab20c",
104+
clim=units_limit
105+
)
106+
pyvista_plotter.add_mesh(
107+
collars,
108+
point_size=10,
109+
render_points_as_spheres=True
110+
)
111+
pyvista_plotter.add_point_labels(
112+
points=borehole_set.collars.collar_loc.points,
113+
labels=borehole_set.collars.ids,
114+
point_size=3,
115+
shape_opacity=0.5,
116+
font_size=12,
117+
bold=True
118+
)
119+
pyvista_plotter.add_actor(gempy_plot.surface_points_actor)
120+
121+
if ADD_GRID := False:
122+
grid_points = pyvista.PointSet(grid.values)
123+
pyvista_plotter.add_mesh(grid_points, color="red", point_size=1)
124+
125+
pyvista_plotter.show()
126+
127+
128+
def _read_spremberg_borehole_set() -> ss.core.geological_formats.BoreholeSet:
129+
from subsurface.core.reader_helpers.readers_data import GenericReaderFilesHelper
130+
from subsurface.core.geological_formats import Survey
131+
from subsurface.core.geological_formats import Collars
132+
from subsurface.core.geological_formats.boreholes.boreholes import MergeOptions
133+
from subsurface.core.geological_formats import BoreholeSet
134+
from subsurface.modules.reader.wells.read_borehole_interface import read_lith, read_survey, read_collar
135+
136+
reader: GenericReaderFilesHelper = GenericReaderFilesHelper(
137+
file_or_buffer=os.getenv("PATH_TO_SPREMBERG_STRATIGRAPHY"),
138+
columns_map={
139+
'hole_id' : 'id',
140+
'depth_from': 'top',
141+
'depth_to' : 'base',
142+
'lit_code' : 'component lith'
143+
}
144+
)
145+
lith: pd.DataFrame = read_lith(reader)
146+
reader: GenericReaderFilesHelper = GenericReaderFilesHelper(
147+
file_or_buffer=os.getenv("PATH_TO_SPREMBERG_SURVEY"),
148+
columns_map={
149+
'depth' : 'md',
150+
'dip' : 'dip',
151+
'azimuth': 'azi'
152+
},
153+
)
154+
df = read_survey(reader)
155+
survey: Survey = Survey.from_df(df, None)
156+
survey.update_survey_with_lith(lith)
157+
reader_collar: GenericReaderFilesHelper = GenericReaderFilesHelper(
158+
file_or_buffer=os.getenv("PATH_TO_SPREMBERG_COLLAR"),
159+
header=0,
160+
usecols=[0, 1, 2, 4],
161+
columns_map={
162+
"hole_id" : "id", # ? Index name is not mapped
163+
"X_GK5_incl_inserted": "x",
164+
"Y__incl_inserted" : "y",
165+
"Z_GK" : "z"
166+
}
167+
)
168+
169+
df_collar = read_collar(reader_collar)
170+
collar = Collars.from_df(df_collar)
171+
borehole_set = BoreholeSet(
172+
collars=collar,
173+
survey=survey,
174+
merge_option=MergeOptions.INTERSECT
175+
)
176+
177+
return borehole_set

0 commit comments

Comments
 (0)