Add example files for gravity and petrophysics inversion (#9)

# Add vector-geology examples and update development logs

This PR adds two complete examples for vector-geology: gravity inversion and petrophysics inversion. These examples demonstrate how to use GemPy for probabilistic modeling in geological contexts.

## Changes:

- Added `.env.example` file with template for API tokens and folder paths
- Updated `.gitignore` to exclude `.nc` files and version file
- Updated development logs with more detailed task tracking for Vector examples
- Added complete example code for gravity inversion with Bayesian methods
- Added complete example code for petrophysics inversion using borehole data
- Added auxiliary functions for both examples in separate files
- Renamed model function in `model_examples.py` to be more descriptive
- Updated test to use the renamed model function
- Added subsurface-terra to dev requirements

The examples include visualization, data processing, and probabilistic modeling techniques that showcase GemPy's capabilities for geological modeling with uncertainty quantification.

Author: Leguark

.env.example

@@ -0,0 +1,2 @@
+LIQUID_EARTH_API_TOKEN = "foobar-1499cbe49808c8cfef14303d2aaacXDx-w1vc1bD065Wj0ABi3HD2fzQ"
+BOREHOLES_SOUTH_FOLDER = "/mnt/d/OneDrive - .../boreholes_south"

.gitignore

@@ -2,6 +2,7 @@
 __pycache__/
 *.py[cod]
 *$py.class
+*.nc
 
 # C extensions
 *.so
@@ -160,4 +161,5 @@ cython_debug/
 .idea/
 
 docs/build
-docs/source
+docs/source
+/gempy_probability/_version.py

docs/dev_logs/2025_May.md

@@ -5,12 +5,17 @@
   - [ ] Keep extracting code into modules and functions 
   - [ ] Posterior analysis for gempy 
 - [ ] Investigate speed
-- [ ] Investigate segmentation function
+- [x] Investigate segmentation function
   https://chatgpt.com/share/68172f53-4f7c-8000-8f42-205a033ed899
 - [ ] Create a robust api to make sure we do not mess up assigning priors to gempy parameters
 - [ ] do we need to use index_put or just normal indexing?
 --- Vector examples
-- [ ] Bring the Vector example here...
+- [ ] Bring the Vector examples here...
+  - Stonepark
+    - Gravity
+    - Multidimensional likelihood
+  - spremberg
+    - Categories
 - [ ] Using categories as likelihood function (maybe I need to check out with Sam)
 --- MinEye
 - [ ] Revive magnetics
@@ -20,11 +25,13 @@
 - [x] Compressing code
 - [x] GemPy posterior visualization
 - [x] Add proper thickness likelihood
-- [ ] Speed
+- [x] Segmentation function 
+- [ ] @Look-Bellow What would I want to have to do StonePark example as flowless as possible? 
+- [ ] Speed: Blocked by having a model that is not too big.
   - [ ] Time test for torch CPU
   - [ ] Time test for keops CPU
   - [ ] Time test for keops GPU
-  - [ ] Options and config presets
+  - [ ] Options and config presets                                                  
 
 ### Possible Design:
 - Geological Model
@@ -38,4 +45,13 @@
 
 ### Likelihood functions
 - ? Should I start to create a catalog of likelihood functions that goes with gempy?
-  - This would be useful to be able to visualize whatever the likelihood function represents independent of the probabilistic model.
+  - This would be useful to be able to visualize whatever the likelihood function represents independent of the probabilistic model.
+  - 
+### What would I want to have to do StonePark example as flowless as possible? 
+1. A way to save the model already with the correct **nugget effects** instead of having to "load nuggets" (TODO in gempy)
+2. Likelihood function plug and play
+3. Some sort of pre-analysis if the model looks good. Maybe some metrics about how long it takes to run the model. 
+   1. Condition number
+   2. Acceptance rate of the posterior
+4. Easy pykeops setting using .env
+5. Some sort of easy set-up for having a reduced number of priors

examples/examples/README.rst

@@ -0,0 +1,6 @@
+Examples
+========
+
+This examples intend to be a gallery of what other people has done with GemPy.
+Hopefully with the help of the community we can create a complete library of
+structures that can serve as template for new users.

examples/examples/__init__.py

@@ -0,0 +1,119 @@
+
+import numpy as np
+
+import gempy as gp
+import xarray as xr
+from vector_geology.model_building_functions import optimize_nuggets_for_group
+
+
+
+
+def initialize_geo_model(structural_elements: list[gp.data.StructuralElement], extent: list[float],
+                         topography: xr.Dataset, load_nuggets: bool = False
+                         ) -> gp.data.GeoModel:
+    structural_group_red = gp.data.StructuralGroup(
+        name="Red",
+        elements=[structural_elements[i] for i in [0, 4, 8]],
+        structural_relation=gp.data.StackRelationType.ERODE
+    )
+
+    # Any, Probably we can decimize this an extra notch
+    structural_group_green = gp.data.StructuralGroup(
+        name="Green",
+        elements=[structural_elements[i] for i in [5]],
+        structural_relation=gp.data.StackRelationType.ERODE
+    )
+
+    # Blue range 2 cov 4
+    structural_group_blue = gp.data.StructuralGroup(
+        name="Blue",
+        elements=[structural_elements[i] for i in [2, 3]],
+        structural_relation=gp.data.StackRelationType.ERODE
+    )
+
+    structural_group_intrusion = gp.data.StructuralGroup(
+        name="Intrusion",
+        elements=[structural_elements[i] for i in [1]],
+        structural_relation=gp.data.StackRelationType.ERODE
+    )
+
+    structural_groups = [structural_group_intrusion, structural_group_green, structural_group_blue, structural_group_red]
+    structural_frame = gp.data.StructuralFrame(
+        structural_groups=structural_groups[2:],
+        color_gen=gp.data.ColorsGenerator()
+    )
+    # TODO: If elements do not have color maybe loop them on structural frame constructor?
+
+    geo_model: gp.data.GeoModel = gp.create_geomodel(
+        project_name='Tutorial_ch1_1_Basics',
+        extent=extent,
+        resolution=[20, 10, 20],
+        refinement=5,  # * Here we define the number of octree levels. If octree levels are defined, the resolution is ignored.
+        structural_frame=structural_frame
+    )
+
+    if topography is not None:
+        gp.set_topography_from_arrays(
+            grid=geo_model.grid,
+            xyz_vertices=topography.vertex.values
+        )
+
+    if load_nuggets:
+        import os
+
+        project_root = os.getcwd()
+        path_to_temp = os.path.join(project_root, "../temp")
+        apply_optimized_nuggets(
+            geo_model=geo_model,
+            loaded_nuggets_red=(np.load(path_to_temp + "/nuggets_Red.npy")),
+            loaded_nuggets_blue=(np.load(path_to_temp + "/nuggets_Blue.npy")),
+            loaded_nuggets_green=(np.load(path_to_temp + "/nuggets_Green.npy"))
+        )
+
+    geo_model.structural_frame.get_element_by_name("KKR").color = "#A46283"
+    geo_model.structural_frame.get_element_by_name("LGR").color = "#6394A4"
+    geo_model.structural_frame.get_element_by_name("WAL").color = "#72A473"
+    geo_model.structural_frame.get_element_by_name("ABL").color = "#1D3943"
+    geo_model.structural_frame.basement_color = "#8B4220"
+
+    geo_model.update_transform()
+
+    return geo_model
+
+def setup_geophysics(env_path: str, geo_model: gp.data.GeoModel):
+    import pandas as pd
+    from dotenv import dotenv_values
+    config = dotenv_values()
+
+    df = pd.read_csv(
+        filepath_or_buffer=config.get(env_path),
+        sep=',',
+        header=0
+    )
+
+    # Remove the items that have X > 5650000
+    df = df[df['X'] < 565000]
+
+    interesting_columns = df[['X', 'Y', 'Bouguer_267_complete']]
+    # %%
+    device_location = interesting_columns[['X', 'Y']]
+    # stack 0 to the z axis
+    device_location['Z'] = 0
+
+    gp.set_centered_grid(
+        grid=geo_model.grid,
+        centers=device_location,
+        resolution=np.array([10, 10, 15]),
+        radius=np.array([5000, 5000, 5000])
+    )
+
+    gravity_gradient = gp.calculate_gravity_gradient(geo_model.grid.centered_grid)
+
+    # %%
+    from gempy_engine.core.backend_tensor import BackendTensor
+    geo_model.geophysics_input = gp.data.GeophysicsInput(
+        tz=BackendTensor.t.array(gravity_gradient),
+        densities=BackendTensor.t.array([2.61, 2.92, 3.1, 2.92, 2.61, 2.61]),
+    )
+
+    return interesting_columns

examples/examples/_aux_func_gravity_inversion.py

@@ -0,0 +1,94 @@
+import os
+import xarray as xr
+import numpy as np
+import gempy as gp
+import subsurface
+
+
+def process_file(filename, global_extent, color_generator: gp.data.ColorsGenerator):
+    dataset: xr.Dataset = xr.open_dataset(filename)
+    file_extent = _calculate_extent(dataset)
+    global_extent = _update_global_extent(global_extent, file_extent)
+
+    base, ext = os.path.splitext(filename)
+    base = os.path.basename(base)
+    structural_element = _extract_surface_points_and_orientations(dataset, base, color_generator)
+
+    return structural_element, global_extent
+
+
+# Function to calculate the extent (ROI) for a single dataset
+def _calculate_extent(dataset):
+    x_coord, y_coord, z_coord = dataset.vertex[:, 0], dataset.vertex[:, 1], dataset.vertex[:, 2]
+    return [x_coord.min().values, x_coord.max().values, y_coord.min().values, y_coord.max().values, z_coord.min().values, z_coord.max().values]
+
+
+# Function to extract surface points and orientations from the dataset
+
+
+def _extract_surface_points_and_orientations(dataset, name, color_generator) -> gp.data.StructuralElement:
+    name = name.replace("Stonepark_", "")
+
+    # Extract surface points and orientations
+    unstruct = subsurface.UnstructuredData(dataset)
+    ts = subsurface.TriSurf(mesh=unstruct)
+    triangulated_mesh = subsurface.visualization.to_pyvista_mesh(ts)
+    # print(name)
+    # subsurface.visualization.pv_plot([triangulated_mesh])
+
+    # Decimate the mesh and compute normals
+    decimated_mesh = triangulated_mesh.decimate_pro(0.95)
+    normals = decimated_mesh.compute_normals(point_normals=False, cell_normals=True, consistent_normals=True)
+    normals_array = np.array(normals.cell_data["Normals"])
+
+    # Filter points
+    vertex_sp = normals.cell_centers().points
+    vertex_grads = normals_array
+    positive_z = vertex_grads[:, 2] > 0
+    largest_component_z = np.all(vertex_grads[:, 2:3] >= np.abs(vertex_grads[:, :2]), axis=1)
+    filter_mask = np.logical_and(positive_z, largest_component_z)
+
+    surface_points_xyz = vertex_sp[filter_mask]
+    orientations_gxyz = vertex_grads[filter_mask]
+    nuggets = np.ones(len(surface_points_xyz)) * 0.000001
+
+    # Create SurfacePointsTable and OrientationsTable
+    every = 10
+    surface_points = gp.data.SurfacePointsTable.from_arrays(
+        x=surface_points_xyz[::every, 0],
+        y=surface_points_xyz[::every, 1],
+        z=surface_points_xyz[::every, 2],
+        names=name,
+        nugget=nuggets[::every]
+    )
+
+    every = 10
+    orientations = gp.data.OrientationsTable.from_arrays(
+        x=surface_points_xyz[::every, 0],
+        y=surface_points_xyz[::every, 1],
+        z=surface_points_xyz[::every, 2],
+        G_x=orientations_gxyz[::every, 0],
+        G_y=orientations_gxyz[::every, 1],
+        G_z=orientations_gxyz[::every, 2],
+        nugget=4,
+        names=name
+    )
+
+    structural_element = gp.data.StructuralElement(
+        name=name,
+        surface_points=surface_points,
+        orientations=orientations,
+        color=next(color_generator)
+    )
+
+    return structural_element
+
+
+# Function to update the global extent based on each file's extent
+def _update_global_extent(global_extent, file_extent):
+    if not global_extent:
+        return file_extent
+    else:
+        return [min(global_extent[0], file_extent[0]), max(global_extent[1], file_extent[1]),
+                min(global_extent[2], file_extent[2]), max(global_extent[3], file_extent[3]),
+                min(global_extent[4], file_extent[4]), max(global_extent[5], file_extent[5])]