.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_basic/03_gravity_forward_model.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_basic_03_gravity_forward_model.py: Forward Gravity Modeling - Tharsis Region ========================================== This example demonstrates forward gravity modeling using a 3D geological model. We compute the gravity response of the geological structure at real measurement locations. **Overview** Forward gravity modeling predicts the gravitational field anomaly caused by subsurface density contrasts. This is a crucial tool in geophysical exploration for: * Validating geological interpretations against field measurements * Identifying density anomalies associated with ore deposits or intrusions * Constraining geological models before inversion * Understanding the relationship between geology and geophysical signatures **Geophysical Background** Gravity surveys measure tiny variations in Earth's gravitational field caused by lateral density variations in the subsurface. The gravity anomaly (Δg) at a point depends on: * The **density contrast** between geological units (Δρ) * The **volume** and **geometry** of anomalous bodies * The **distance** from the observation point to the density anomaly Typical units: * mGal (milligal) = 10⁻⁵ m/s² = 10⁻³ cm/s² * μGal (microgal) = 10⁻⁸ m/s² = 10⁻⁶ cm/s² Modern gravimeters can measure to ~10 μGal precision. **The Tharsis Case Study** The Tharsis plutonite intrusion (density ~2.9 g/cm³) intruded into sedimentary host rocks (density ~2.3 g/cm³), creating a positive gravity anomaly. This density contrast of ~0.6 g/cm³ produces a measurable signal that we can model and compare with observations. **Workflow Summary** 1. Create geological model (similar to Example 01) 2. Load field gravity measurements 3. Set up computation grids at measurement locations 4. Assign density values to geological units 5. Compute forward gravity response 6. Compare with observations and analyze residuals .. GENERATED FROM PYTHON SOURCE LINES 49-53 .. code-block:: Python import dotenv from mineye.GeoModel.helper_plotter import plot_model_and_gravity_sensors .. GENERATED FROM PYTHON SOURCE LINES 54-72 Import Libraries ---------------- **Core Libraries**: * **NumPy**: Numerical array operations * **PyTorch**: Tensor operations with automatic differentiation (useful for inversions) * **Matplotlib**: Plotting and visualization **Geological Modeling**: * **GemPy**: 3D implicit geological modeling framework * **GeoPandas**: Geospatial data handling (for loading gravity measurements with coordinates) **Mineye Utilities**: * ``align_forward_to_observed``: Aligns modeled gravity to observed data statistics * ``normalize``: Data normalization and preprocessing functions .. GENERATED FROM PYTHON SOURCE LINES 72-88 .. code-block:: Python dotenv.load_dotenv() import numpy as np import torch from matplotlib import pyplot as plt import gempy as gp import geopandas as gpd from mineye.GeoModel.geophysics import align_forward_to_observed from mineye.GeoModel.model_one.probabilistic_model import normalize # Set random seed for reproducibility np.random.seed(1234) .. GENERATED FROM PYTHON SOURCE LINES 89-90 Import paths configuration .. GENERATED FROM PYTHON SOURCE LINES 90-92 .. code-block:: Python from mineye.config import paths .. GENERATED FROM PYTHON SOURCE LINES 93-96 Define Model Extent and Resolution ----------------------------------- The model covers the Tharsis mining district .. GENERATED FROM PYTHON SOURCE LINES 96-109 .. code-block:: Python min_x = -707_521 # * Cropping the corrupted area of the geotiff max_x = -675558 min_y = 4526832 max_y = 4551949 max_z = 505 model_depth = -500 extent = [min_x, max_x, min_y, max_y, model_depth, max_z] # Model resolution: use octree with refinement level 5 resolution = None # Using octrees instead of regular grid refinement = 5 .. GENERATED FROM PYTHON SOURCE LINES 110-113 Get Data Paths -------------- Load paths to structural and geophysical data .. GENERATED FROM PYTHON SOURCE LINES 113-122 .. code-block:: Python mod_or_path = paths.get_orientations_path() mod_pts_path = paths.get_points_path() gravity_data_path = paths.get_gravity_data_path() print(f"Orientations: {mod_or_path}") print(f"Points: {mod_pts_path}") print(f"Gravity data: {gravity_data_path}") .. rst-class:: sphx-glr-script-out .. code-block:: none Orientations: /home/leguark/PycharmProjects/Mineye-Terranigma/examples/Data/Model_Input_Data/Simple-Models/orientations_mod.csv Points: /home/leguark/PycharmProjects/Mineye-Terranigma/examples/Data/Model_Input_Data/Simple-Models/points_mod.csv Gravity data: /home/leguark/PycharmProjects/Mineye-Terranigma/examples/Data/General_Input_Data/Geophysical_Cleaned_Data/cleaned_gravity_data.geojson .. GENERATED FROM PYTHON SOURCE LINES 123-126 Create GemPy Geological Model ------------------------------ Build the model structure with imported structural data .. GENERATED FROM PYTHON SOURCE LINES 126-138 .. code-block:: Python simple_geo_model = gp.create_geomodel( project_name='gravity_model', extent=extent, refinement=refinement, resolution=resolution, importer_helper=gp.data.ImporterHelper( path_to_orientations=mod_or_path, path_to_surface_points=mod_pts_path, ) ) .. GENERATED FROM PYTHON SOURCE LINES 139-142 Map Geological Units -------------------- Define the stratigraphic stack .. GENERATED FROM PYTHON SOURCE LINES 142-150 .. code-block:: Python gp.map_stack_to_surfaces( gempy_model=simple_geo_model, mapping_object={ "Tournaisian_Plutonites": ["Tournaisian Plutonites"], } ) .. raw:: html
Structural Groups: StructuralGroup:
Name:Tournaisian_Plutonites
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Tournaisian Plutonites
Fault Relations:
Tournaisia...
Tournaisian_Plutonites
True
False


.. GENERATED FROM PYTHON SOURCE LINES 151-174 Load Gravity Observations -------------------------- Read actual gravity measurements from the field **Gravity Survey Data**: Gravity surveys involve systematic measurements of gravitational acceleration across an area of interest. The data typically includes: * **Spatial coordinates** (X, Y, Z): Location of each measurement station * **Bouguer anomaly**: Gravity after correcting for latitude, elevation, terrain, and tidal effects * **Free-air anomaly**: Simpler correction (elevation only) * **Measurement metadata**: Date, instrument ID, drift corrections, etc. The Bouguer anomaly (``VALU_BOU267``) represents the gravity signal attributed purely to subsurface density variations, making it ideal for geological interpretation. **Data Quality Considerations**: * Survey spacing and coverage affect resolution * Measurement precision (typically 10-50 μGal for ground surveys) * Systematic errors from incomplete corrections * Regional trends vs. local anomalies .. GENERATED FROM PYTHON SOURCE LINES 174-181 .. code-block:: Python gravity_data = gpd.read_file(gravity_data_path) observed_gravity = gravity_data['VALU_BOU267'].values # in mGal print(f"Number of gravity observations: {len(observed_gravity)}") print(f"Gravity range: {observed_gravity.min():.2f} to {observed_gravity.max():.2f} mGal") .. rst-class:: sphx-glr-script-out .. code-block:: none Number of gravity observations: 260 Gravity range: 9.51 to 47.32 mGal .. GENERATED FROM PYTHON SOURCE LINES 182-197 Set Up Gravity Measurement Locations ------------------------------------- Create measurement grid at actual device locations **Centered Grid Approach**: For efficient gravity computation, GemPy uses "centered grids" - small 3D volumes around each measurement point. This approach: * Focuses computation only where observations exist (not the entire model volume) * Maintains high resolution near observation points * Significantly reduces computational cost for large numbers of stations Each measurement location becomes the center of a local grid where the density distribution is evaluated in detail. .. GENERATED FROM PYTHON SOURCE LINES 197-206 .. code-block:: Python xy_ravel = np.column_stack([ np.array(gravity_data.geometry.x.values), np.array(gravity_data.geometry.y.values), np.full(len(gravity_data), 500) # Set Z to surface elevation ]) print(f"Using {len(xy_ravel)} actual measurement points") .. rst-class:: sphx-glr-script-out .. code-block:: none Using 260 actual measurement points .. GENERATED FROM PYTHON SOURCE LINES 207-246 Configure Centered Grid for Gravity Computation ------------------------------------------------ Set up centered grid around each measurement point **Grid Parameters**: * **centers**: Coordinates of measurement stations * **resolution**: Number of cells in each direction [X, Y, Z] - [10, 10, 15] creates 10×10×15 = 1,500 cells per station * **radius**: Extent of local grid in meters [X, Y, Z] - [5000, 5000, 5000] = 5 km radius in all directions The grid resolution balances accuracy vs. computation time. Higher resolution captures finer details of the density distribution but increases computation cost. **The 1/r² Nature of Gravity** Gravitational attraction follows an inverse-square law: .. math:: F = G \frac{m_1 m_2}{r^2} This means that the contribution of a mass element to the gravity measured at a station decreases rapidly with distance. **Why the Radius Matters in set_centered_grid** The `radius` parameter in `set_centered_grid` defines the 3D volume around each station where GemPy will evaluate the density distribution. * **Too small radius**: Fails to capture the signal from larger or deeper structures, leading to an underestimated or "chopped" gravity response. * **Too large radius**: Increases computational cost without significantly improving accuracy, as distant cells have negligible influence. In this example, a 5 km radius ([5000, 5000, 5000]) is chosen to encompass the entire Tharsis plutonite intrusion, ensuring its full gravitational signature is captured. .. GENERATED FROM PYTHON SOURCE LINES 246-254 .. code-block:: Python gp.set_centered_grid( grid=simple_geo_model.grid, centers=xy_ravel, resolution=np.array([10, 10, 15]), radius=np.array([5000, 5000, 5000]) ) .. rst-class:: sphx-glr-script-out .. code-block:: none Active grids: GridTypes.OCTREE|CENTERED|NONE CenteredGrid(centers=array([[-7.04570250e+05, 4.54824444e+06, 5.00000000e+02], [-7.04382532e+05, 4.53206604e+06, 5.00000000e+02], [-6.90136472e+05, 4.54443298e+06, 5.00000000e+02], [-7.07002148e+05, 4.53909364e+06, 5.00000000e+02], [-7.04699537e+05, 4.53196859e+06, 5.00000000e+02], [-7.00950961e+05, 4.54679353e+06, 5.00000000e+02], [-6.94189651e+05, 4.54889852e+06, 5.00000000e+02], [-7.04846430e+05, 4.53911004e+06, 5.00000000e+02], [-7.03528058e+05, 4.53677639e+06, 5.00000000e+02], [-6.97783816e+05, 4.53369043e+06, 5.00000000e+02], [-7.00701909e+05, 4.53205343e+06, 5.00000000e+02], [-6.89387223e+05, 4.54888667e+06, 5.00000000e+02], [-7.06356826e+05, 4.54338128e+06, 5.00000000e+02], [-6.99725791e+05, 4.53900142e+06, 5.00000000e+02], [-7.04880761e+05, 4.53972244e+06, 5.00000000e+02], [-6.93539735e+05, 4.53356575e+06, 5.00000000e+02], [-6.89703778e+05, 4.54411549e+06, 5.00000000e+02], [-6.95657885e+05, 4.54849206e+06, 5.00000000e+02], [-7.01340437e+05, 4.54453812e+06, 5.00000000e+02], [-7.03113495e+05, 4.54434831e+06, 5.00000000e+02], [-6.96260288e+05, 4.54253904e+06, 5.00000000e+02], [-7.05269751e+05, 4.53986829e+06, 5.00000000e+02], [-7.02942617e+05, 4.53742719e+06, 5.00000000e+02], [-7.01668401e+05, 4.53809354e+06, 5.00000000e+02], [-6.96778079e+05, 4.53456324e+06, 5.00000000e+02], [-6.90189586e+05, 4.53590404e+06, 5.00000000e+02], [-6.94278014e+05, 4.54039014e+06, 5.00000000e+02], [-7.05932225e+05, 4.54103010e+06, 5.00000000e+02], [-6.94347590e+05, 4.53669093e+06, 5.00000000e+02], [-6.94530807e+05, 4.53812585e+06, 5.00000000e+02], [-7.06304300e+05, 4.53856975e+06, 5.00000000e+02], [-6.99847474e+05, 4.53020227e+06, 5.00000000e+02], [-7.00873497e+05, 4.54048111e+06, 5.00000000e+02], [-6.95459414e+05, 4.53584489e+06, 5.00000000e+02], [-7.07148489e+05, 4.54089640e+06, 5.00000000e+02], [-6.97998064e+05, 4.53419107e+06, 5.00000000e+02], [-6.91595188e+05, 4.53870027e+06, 5.00000000e+02], [-7.06430337e+05, 4.53622273e+06, 5.00000000e+02], [-6.99279789e+05, 4.54864673e+06, 5.00000000e+02], [-7.00033658e+05, 4.52899618e+06, 5.00000000e+02], [-7.05136851e+05, 4.53298811e+06, 5.00000000e+02], [-7.01343701e+05, 4.53146698e+06, 5.00000000e+02], [-6.94948029e+05, 4.53785148e+06, 5.00000000e+02], [-7.02561885e+05, 4.54905334e+06, 5.00000000e+02], [-6.93552418e+05, 4.54198266e+06, 5.00000000e+02], [-6.99527196e+05, 4.54773581e+06, 5.00000000e+02], [-7.01848423e+05, 4.54422076e+06, 5.00000000e+02], [-7.01110480e+05, 4.54442829e+06, 5.00000000e+02], [-6.97729207e+05, 4.54472250e+06, 5.00000000e+02], [-6.97146948e+05, 4.54465158e+06, 5.00000000e+02], [-7.02279941e+05, 4.54567818e+06, 5.00000000e+02], [-7.02500434e+05, 4.54619937e+06, 5.00000000e+02], [-7.01558235e+05, 4.54584715e+06, 5.00000000e+02], [-7.02871591e+05, 4.54468798e+06, 5.00000000e+02], [-7.03700811e+05, 4.54459757e+06, 5.00000000e+02], [-7.04147383e+05, 4.54440278e+06, 5.00000000e+02], [-7.07131689e+05, 4.54439096e+06, 5.00000000e+02], [-7.06458744e+05, 4.54705891e+06, 5.00000000e+02], [-7.05366723e+05, 4.54846972e+06, 5.00000000e+02], [-7.05436911e+05, 4.54806242e+06, 5.00000000e+02], [-7.05547079e+05, 4.54714491e+06, 5.00000000e+02], [-6.92878521e+05, 4.54534566e+06, 5.00000000e+02], [-6.93672250e+05, 4.54458621e+06, 5.00000000e+02], [-6.94720653e+05, 4.54491513e+06, 5.00000000e+02], [-6.90842088e+05, 4.54549207e+06, 5.00000000e+02], [-6.91144549e+05, 4.54609948e+06, 5.00000000e+02], [-6.88911827e+05, 4.54541877e+06, 5.00000000e+02], [-6.92466066e+05, 4.54912199e+06, 5.00000000e+02], [-6.92675208e+05, 4.54880892e+06, 5.00000000e+02], [-6.94114991e+05, 4.54838209e+06, 5.00000000e+02], [-6.90032797e+05, 4.54769453e+06, 5.00000000e+02], [-6.90656637e+05, 4.54794335e+06, 5.00000000e+02], [-6.95358046e+05, 4.54895682e+06, 5.00000000e+02], [-6.93118147e+05, 4.54691561e+06, 5.00000000e+02], [-6.91034221e+05, 4.54841763e+06, 5.00000000e+02], [-7.01197862e+05, 4.54795971e+06, 5.00000000e+02], [-6.95771047e+05, 4.54708854e+06, 5.00000000e+02], [-6.95592360e+05, 4.54652510e+06, 5.00000000e+02], [-6.94909127e+05, 4.54682550e+06, 5.00000000e+02], [-6.94523994e+05, 4.54804049e+06, 5.00000000e+02], [-6.92347602e+05, 4.54644112e+06, 5.00000000e+02], [-6.98635862e+05, 4.54619901e+06, 5.00000000e+02], [-6.99251249e+05, 4.54490114e+06, 5.00000000e+02], [-7.00027786e+05, 4.54698477e+06, 5.00000000e+02], [-6.99853919e+05, 4.54746797e+06, 5.00000000e+02], [-6.98898971e+05, 4.54787341e+06, 5.00000000e+02], [-6.98433212e+05, 4.54748884e+06, 5.00000000e+02], [-6.98246166e+05, 4.54712493e+06, 5.00000000e+02], [-6.98987687e+05, 4.54466052e+06, 5.00000000e+02], [-6.98034309e+05, 4.54481454e+06, 5.00000000e+02], [-7.03419465e+05, 4.54643054e+06, 5.00000000e+02], [-7.03097299e+05, 4.54700280e+06, 5.00000000e+02], [-7.02423402e+05, 4.54795682e+06, 5.00000000e+02], [-7.00026697e+05, 4.54807220e+06, 5.00000000e+02], [-7.03474125e+05, 4.53288773e+06, 5.00000000e+02], [-7.03810666e+05, 4.53255646e+06, 5.00000000e+02], [-7.04673247e+05, 4.53436305e+06, 5.00000000e+02], [-7.03698423e+05, 4.53427360e+06, 5.00000000e+02], [-7.05245638e+05, 4.53328072e+06, 5.00000000e+02], [-7.02630887e+05, 4.53936393e+06, 5.00000000e+02], [-7.03702556e+05, 4.54003708e+06, 5.00000000e+02], [-7.06361499e+05, 4.54134820e+06, 5.00000000e+02], [-7.04223671e+05, 4.53849332e+06, 5.00000000e+02], [-7.04356379e+05, 4.53767727e+06, 5.00000000e+02], [-7.04665012e+05, 4.53933458e+06, 5.00000000e+02], [-7.06746112e+05, 4.53924229e+06, 5.00000000e+02], [-7.07058647e+05, 4.54248477e+06, 5.00000000e+02], [-7.06766004e+05, 4.54270708e+06, 5.00000000e+02], [-7.03396645e+05, 4.54135021e+06, 5.00000000e+02], [-7.06495218e+05, 4.54237312e+06, 5.00000000e+02], [-7.02357301e+05, 4.54027131e+06, 5.00000000e+02], [-7.01650512e+05, 4.54085429e+06, 5.00000000e+02], [-7.02891561e+05, 4.54357107e+06, 5.00000000e+02], [-7.07191497e+05, 4.54132345e+06, 5.00000000e+02], [-7.06988196e+05, 4.53972932e+06, 5.00000000e+02], [-7.02039309e+05, 4.54268183e+06, 5.00000000e+02], [-7.03405547e+05, 4.53566539e+06, 5.00000000e+02], [-7.03266909e+05, 4.53531292e+06, 5.00000000e+02], [-7.00236264e+05, 4.54215213e+06, 5.00000000e+02], [-7.03098064e+05, 4.53675239e+06, 5.00000000e+02], [-7.00548764e+05, 4.54083495e+06, 5.00000000e+02], [-7.01069088e+05, 4.53817433e+06, 5.00000000e+02], [-7.02545850e+05, 4.53694946e+06, 5.00000000e+02], [-7.01127016e+05, 4.53953768e+06, 5.00000000e+02], [-7.01413474e+05, 4.53795171e+06, 5.00000000e+02], [-7.01057167e+05, 4.53751050e+06, 5.00000000e+02], [-7.00774764e+05, 4.53420618e+06, 5.00000000e+02], [-6.99526050e+05, 4.53657202e+06, 5.00000000e+02], [-6.98705002e+05, 4.53655399e+06, 5.00000000e+02], [-6.98738611e+05, 4.53596090e+06, 5.00000000e+02], [-6.98138891e+05, 4.53541035e+06, 5.00000000e+02], [-7.00814354e+05, 4.53274251e+06, 5.00000000e+02], [-7.01536976e+05, 4.53316425e+06, 5.00000000e+02], [-6.95709215e+05, 4.53386932e+06, 5.00000000e+02], [-6.96976393e+05, 4.53689850e+06, 5.00000000e+02], [-6.99658441e+05, 4.53302653e+06, 5.00000000e+02], [-6.96528748e+05, 4.53508901e+06, 5.00000000e+02], [-6.96655907e+05, 4.53413284e+06, 5.00000000e+02], [-6.95936549e+05, 4.53555451e+06, 5.00000000e+02], [-6.94863559e+05, 4.53316759e+06, 5.00000000e+02], [-6.95976080e+05, 4.53600434e+06, 5.00000000e+02], [-6.96879191e+05, 4.53336600e+06, 5.00000000e+02], [-6.96445459e+05, 4.53297440e+06, 5.00000000e+02], [-6.97873969e+05, 4.53319816e+06, 5.00000000e+02], [-6.94769803e+05, 4.53381582e+06, 5.00000000e+02], [-6.94221245e+05, 4.53819111e+06, 5.00000000e+02], [-6.93857000e+05, 4.53864873e+06, 5.00000000e+02], [-6.93360409e+05, 4.53400280e+06, 5.00000000e+02], [-6.93017493e+05, 4.53404887e+06, 5.00000000e+02], [-6.93705289e+05, 4.53702628e+06, 5.00000000e+02], [-6.92957109e+05, 4.53674162e+06, 5.00000000e+02], [-6.91953387e+05, 4.53651263e+06, 5.00000000e+02], [-6.94163748e+05, 4.53960164e+06, 5.00000000e+02], [-6.93028813e+05, 4.54287968e+06, 5.00000000e+02], [-6.93133427e+05, 4.54346203e+06, 5.00000000e+02], [-6.97132629e+05, 4.54212935e+06, 5.00000000e+02], [-6.97180785e+05, 4.54260705e+06, 5.00000000e+02], [-6.97719242e+05, 4.54283803e+06, 5.00000000e+02], [-6.90694925e+05, 4.54158545e+06, 5.00000000e+02], [-6.90350442e+05, 4.54163026e+06, 5.00000000e+02], [-6.89598979e+05, 4.54129335e+06, 5.00000000e+02], [-6.89177714e+05, 4.54071305e+06, 5.00000000e+02], [-6.92605683e+05, 4.53506806e+06, 5.00000000e+02], [-6.92447826e+05, 4.53535993e+06, 5.00000000e+02], [-6.89925868e+05, 4.53980363e+06, 5.00000000e+02], [-6.91488423e+05, 4.53483537e+06, 5.00000000e+02], [-6.91400717e+05, 4.53929017e+06, 5.00000000e+02], [-6.91219303e+05, 4.53841402e+06, 5.00000000e+02], [-6.90900320e+05, 4.53735252e+06, 5.00000000e+02], [-6.91760824e+05, 4.53889234e+06, 5.00000000e+02], [-6.92013579e+05, 4.53827256e+06, 5.00000000e+02], [-6.93278788e+05, 4.53938381e+06, 5.00000000e+02], [-6.91108667e+05, 4.53975920e+06, 5.00000000e+02], [-6.95312407e+05, 4.54273825e+06, 5.00000000e+02], [-6.96085901e+05, 4.54297381e+06, 5.00000000e+02], [-6.96506017e+05, 4.54106574e+06, 5.00000000e+02], [-6.89036650e+05, 4.54373753e+06, 5.00000000e+02], [-6.96502841e+05, 4.54035450e+06, 5.00000000e+02], [-6.91995894e+05, 4.54325403e+06, 5.00000000e+02], [-6.95746810e+05, 4.54109351e+06, 5.00000000e+02], [-6.95624725e+05, 4.54201843e+06, 5.00000000e+02], [-6.92006343e+05, 4.54237855e+06, 5.00000000e+02], [-6.91133139e+05, 4.54262171e+06, 5.00000000e+02], [-6.90716496e+05, 4.54251429e+06, 5.00000000e+02], [-6.95074940e+05, 4.54313969e+06, 5.00000000e+02], [-6.99027896e+05, 4.53951537e+06, 5.00000000e+02], [-7.01022098e+05, 4.54410543e+06, 5.00000000e+02], [-7.00368871e+05, 4.54002472e+06, 5.00000000e+02], [-7.00069026e+05, 4.53993514e+06, 5.00000000e+02], [-6.98734314e+05, 4.53996692e+06, 5.00000000e+02], [-6.97556326e+05, 4.54392817e+06, 5.00000000e+02], [-6.96691841e+05, 4.54402186e+06, 5.00000000e+02], [-7.00005633e+05, 4.54324843e+06, 5.00000000e+02], [-7.05393057e+05, 4.54241772e+06, 5.00000000e+02], [-7.03824349e+05, 4.54416189e+06, 5.00000000e+02], [-6.89313757e+05, 4.53575787e+06, 5.00000000e+02], [-6.89391846e+05, 4.54400568e+06, 5.00000000e+02], [-6.90158384e+05, 4.54355574e+06, 5.00000000e+02], [-6.98887192e+05, 4.53830451e+06, 5.00000000e+02], [-6.99887213e+05, 4.53812986e+06, 5.00000000e+02], [-7.06792156e+05, 4.53263109e+06, 5.00000000e+02], [-6.92886936e+05, 4.54415550e+06, 5.00000000e+02], [-6.92149125e+05, 4.54547872e+06, 5.00000000e+02], [-6.96750137e+05, 4.53696218e+06, 5.00000000e+02], [-6.96729668e+05, 4.54642320e+06, 5.00000000e+02], [-6.80457499e+05, 4.52996092e+06, 5.00000000e+02], [-7.05624398e+05, 4.54814706e+06, 5.00000000e+02], [-7.03076852e+05, 4.53811962e+06, 5.00000000e+02], [-6.92756568e+05, 4.53085230e+06, 5.00000000e+02], [-6.88941539e+05, 4.52909583e+06, 5.00000000e+02], [-7.01749402e+05, 4.53120052e+06, 5.00000000e+02], [-7.05629598e+05, 4.54683619e+06, 5.00000000e+02], [-7.00662443e+05, 4.53723290e+06, 5.00000000e+02], [-6.94478080e+05, 4.54018665e+06, 5.00000000e+02], [-6.89850511e+05, 4.54152314e+06, 5.00000000e+02], [-6.94523994e+05, 4.54804049e+06, 5.00000000e+02], [-7.06608130e+05, 4.54813077e+06, 5.00000000e+02], [-6.91079208e+05, 4.54577696e+06, 5.00000000e+02], [-7.06061105e+05, 4.52881513e+06, 5.00000000e+02], [-7.06724221e+05, 4.53104101e+06, 5.00000000e+02], [-7.04821306e+05, 4.52885057e+06, 5.00000000e+02], [-7.01080481e+05, 4.53021960e+06, 5.00000000e+02], [-6.97759999e+05, 4.53020714e+06, 5.00000000e+02], [-6.97560942e+05, 4.52937808e+06, 5.00000000e+02], [-6.95454344e+05, 4.52900695e+06, 5.00000000e+02], [-6.94517920e+05, 4.52943086e+06, 5.00000000e+02], [-6.93415278e+05, 4.53137302e+06, 5.00000000e+02], [-6.94228488e+05, 4.53198213e+06, 5.00000000e+02], [-6.95617778e+05, 4.53008808e+06, 5.00000000e+02], [-6.94136671e+05, 4.52953148e+06, 5.00000000e+02], [-6.93887658e+05, 4.52999358e+06, 5.00000000e+02], [-6.91133796e+05, 4.52980229e+06, 5.00000000e+02], [-6.94161486e+05, 4.53102414e+06, 5.00000000e+02], [-6.89202793e+05, 4.52921231e+06, 5.00000000e+02], [-6.89300434e+05, 4.53056691e+06, 5.00000000e+02], [-7.01111262e+05, 4.53124209e+06, 5.00000000e+02], [-7.01611930e+05, 4.52951173e+06, 5.00000000e+02], [-6.99811930e+05, 4.53124357e+06, 5.00000000e+02], [-7.00231847e+05, 4.53051682e+06, 5.00000000e+02], [-7.03557165e+05, 4.53159373e+06, 5.00000000e+02], [-7.05055274e+05, 4.53097026e+06, 5.00000000e+02], [-7.06108138e+05, 4.53003321e+06, 5.00000000e+02], [-7.04459358e+05, 4.53002878e+06, 5.00000000e+02], [-7.04977759e+05, 4.53237929e+06, 5.00000000e+02], [-6.95321488e+05, 4.53214056e+06, 5.00000000e+02], [-6.97282739e+05, 4.53132216e+06, 5.00000000e+02], [-6.97979602e+05, 4.53082777e+06, 5.00000000e+02], [-6.92196403e+05, 4.53099830e+06, 5.00000000e+02], [-7.03422717e+05, 4.54244623e+06, 5.00000000e+02], [-7.06538457e+05, 4.52890316e+06, 5.00000000e+02], [-6.99215407e+05, 4.53530158e+06, 5.00000000e+02], [-6.89698506e+05, 4.54704240e+06, 5.00000000e+02], [-6.99284130e+05, 4.53633084e+06, 5.00000000e+02], [-6.96499339e+05, 4.53453115e+06, 5.00000000e+02], [-6.82972390e+05, 4.54842446e+06, 5.00000000e+02], [-6.87802644e+05, 4.54067634e+06, 5.00000000e+02], [-6.93430253e+05, 4.54254364e+06, 5.00000000e+02], [-7.01585545e+05, 4.53598599e+06, 5.00000000e+02], [-6.98010147e+05, 4.53061768e+06, 5.00000000e+02], [-7.00367868e+05, 4.54428712e+06, 5.00000000e+02]]), resolution=array([10, 10, 15]), radius=array([5000, 5000, 5000]), kernel_grid_centers=array([[-5000. , -5000. , -300. ], [-5000. , -5000. , -360. ], [-5000. , -5000. , -383.36972966], ..., [ 5000. , 5000. , -3407.68480754], [ 5000. , 5000. , -4618.11403801], [ 5000. , 5000. , -6300. ]], shape=(1936, 3)), left_voxel_edges=array([[1709.43058496, 1709.43058496, -30. ], [1709.43058496, 1709.43058496, -30. ], [1709.43058496, 1709.43058496, -11.68486483], ..., [1709.43058496, 1709.43058496, -435.56428767], [1709.43058496, 1709.43058496, -605.21461523], [1709.43058496, 1709.43058496, -840.942981 ]], shape=(1936, 3)), right_voxel_edges=array([[1709.43058496, 1709.43058496, -30. ], [1709.43058496, 1709.43058496, -11.68486483], [1709.43058496, 1709.43058496, -16.23606704], ..., [1709.43058496, 1709.43058496, -605.21461523], [1709.43058496, 1709.43058496, -840.942981 ], [1709.43058496, 1709.43058496, -840.942981 ]], shape=(1936, 3))) .. GENERATED FROM PYTHON SOURCE LINES 255-275 Calculate Gravity Gradient --------------------------- Compute the gravity gradient (tz component) for the centered grid **Gravity Tensor Components**: The gravity field can be represented as a gradient tensor with 9 components. For vertical gravity (what gravimeters measure), we primarily use: * **g_z** (or **t_z**): Vertical component of gravitational acceleration * This is what ground-based gravimeters measure The gradient calculation determines the geometric contribution of each grid cell to the total gravity at the observation point, based on: * Distance from cell to observation point * Cell volume * Orientation (vertical direction) These gradients are later multiplied by density values to get the final gravity response. .. GENERATED FROM PYTHON SOURCE LINES 275-279 .. code-block:: Python gravity_gradient = gp.calculate_gravity_gradient(simple_geo_model.grid.centered_grid) print(f"Gravity gradient tensor shape: {gravity_gradient.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Gravity gradient tensor shape: (1936,) .. GENERATED FROM PYTHON SOURCE LINES 280-308 Set Density Values ------------------ Define densities for geological units **Rock Densities**: Density contrasts drive gravity anomalies. Typical rock densities: * **Plutonic rocks** (granite, diorite): 2.60-2.95 g/cm³ * **Sedimentary rocks** (shale, sandstone): 2.20-2.70 g/cm³ * **Volcanic rocks** (basalt): 2.70-3.10 g/cm³ * **Ore deposits** (massive sulfides): 3.50-5.00 g/cm³ **Tharsis Densities**: * **Tournaisian Plutonites**: 2.9 g/cm³ (dense intrusive rocks) * **Devonian Sedimentary Host**: 2.3 g/cm³ (lower density host rocks) * **Density contrast**: +0.6 g/cm³ (creates positive gravity anomaly) These values are based on: * Laboratory measurements of rock samples * Literature values for similar rock types * Calibration with gravity data (in inverse modeling) .. note:: Density uncertainties significantly impact modeling results. Typical uncertainties are ±0.05-0.10 g/cm³, which should be considered in probabilistic inversions. .. GENERATED FROM PYTHON SOURCE LINES 308-321 .. code-block:: Python density_plutonites = 2.9 # g/cm³ density_sedimentary_host = 2.3 # g/cm³ simple_geo_model.geophysics_input = gp.data.GeophysicsInput( tz=gravity_gradient, densities=np.array([density_plutonites, density_sedimentary_host]) ) print(f"Plutonite density: {density_plutonites} g/cm³") print(f"Sedimentary host density: {density_sedimentary_host} g/cm³") print(f"Density contrast: {density_plutonites - density_sedimentary_host:.1f} g/cm³") .. rst-class:: sphx-glr-script-out .. code-block:: none Plutonite density: 2.9 g/cm³ Sedimentary host density: 2.3 g/cm³ Density contrast: 0.6 g/cm³ .. GENERATED FROM PYTHON SOURCE LINES 322-346 Compute Forward Gravity Model ------------------------------ Run the interpolation and gravity computation **Forward Modeling Process**: 1. **Geological interpolation**: Determine which geological unit occupies each grid cell 2. **Density assignment**: Assign the appropriate density to each cell based on lithology 3. **Gravity summation**: Sum the gravitational contribution of all cells using: .. math:: g_z = G \\sum_{i} \\frac{\\Delta\\rho_i \\cdot V_i \\cdot z_i}{r_i^3} where: * G = gravitational constant (6.674×10⁻¹¹ m³/(kg·s²)) * Δρᵢ = density of cell i * Vᵢ = volume of cell i * zᵢ = vertical distance from cell to observation point * rᵢ = total distance from cell to observation point 4. **Unit conversion**: Convert to mGal or μGal The ``mesh_extraction=False`` flag skips 3D mesh generation (not needed for gravity). .. GENERATED FROM PYTHON SOURCE LINES 346-364 .. code-block:: Python simple_geo_model.interpolation_options.mesh_extraction = True gp.set_topography_from_file( grid=simple_geo_model.grid, filepath=paths.get_topography_path(), crop_to_extent=[ simple_geo_model.grid.extent[0], simple_geo_model.grid.extent[2], simple_geo_model.grid.extent[1], simple_geo_model.grid.extent[3] ] ) sol = gp.compute_model(simple_geo_model, validate_serialization=False) print("✓ Forward gravity model computed successfully") .. rst-class:: sphx-glr-script-out .. code-block:: none Active grids: GridTypes.OCTREE|TOPOGRAPHY|CENTERED|NONE Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces ✓ Forward gravity model computed successfully .. GENERATED FROM PYTHON SOURCE LINES 365-384 Extract Gravity Results ----------------------- Get the computed gravity response **Data Alignment**: The ``align_forward_to_observed`` function adjusts the forward model to match the statistical properties of observations. This accounts for: * **Regional trends**: Long-wavelength variations not captured by the local model * **Reference level differences**: Arbitrary datums in gravity data * **Systematic biases**: Uncorrected effects or model simplifications The normalization uses the ``align_to_reference`` method with a 30% extrapolation buffer to handle edge effects. .. note:: This alignment is a practical necessity for forward modeling but should be carefully considered in inverse problems where we're trying to fit the data. .. GENERATED FROM PYTHON SOURCE LINES 384-399 .. code-block:: Python observed_gravity_ugal = observed_gravity * 1000 # Convert mGal to μGal norm_params = normalize( sol.gravity, observed_gravity_ugal, method="align_to_reference", extrapolation_buffer=0.3 ) grav = align_forward_to_observed(sol.gravity, norm_params) print(f"\nComputed gravity values:") print(f" Shape: {grav.shape}") print(f" Range: {grav.min():.2f} to {grav.max():.2f} μGal") print(f" Mean: {grav.mean():.2f} μGal") .. rst-class:: sphx-glr-script-out .. code-block:: none Computing align_to_reference alignment parameters... Alignment params: {'method': 'align_to_reference', 'reference_mean': 25587.184615384616, 'reference_std': 8369.297471927344, 'baseline_forward_mean': -388.1267789929416, 'baseline_forward_std': 1.8054178450074378} Computed gravity values: Shape: torch.Size([260]) Range: -50607.12 to 29982.25 μGal Mean: 25587.18 μGal .. GENERATED FROM PYTHON SOURCE LINES 400-425 Compare with Observations ------------------------- Calculate residuals between observed and computed gravity **Residual Analysis**: Residuals (observed - modeled) reveal how well the geological model explains the data: * **Mean residual**: Should be close to zero after alignment; systematic bias indicates regional trends or model errors * **Standard deviation**: Measures scatter; affected by model accuracy and data noise * **RMSE** (Root Mean Square Error): Overall measure of fit quality .. math:: RMSE = \\sqrt{\\frac{1}{N} \\sum_{i=1}^{N} (g_{obs,i} - g_{mod,i})^2} **Interpreting Residuals**: * **Systematic patterns**: Suggest missing geological features or incorrect geometry * **Random scatter**: Dominated by measurement noise or small-scale geology * **Large residuals**: May indicate: - Incorrect density values - Unmodeled geological structures - Data quality issues - 3D effects not captured by 2D modeling .. GENERATED FROM PYTHON SOURCE LINES 425-436 .. code-block:: Python if isinstance(grav, torch.Tensor): grav = grav.detach().numpy() residuals = observed_gravity_ugal - grav print(f"\nGravity residuals:") print(f" Mean: {residuals.mean():.2f} μGal") print(f" Std: {residuals.std():.2f} μGal") print(f" RMS: {np.sqrt(np.mean(residuals ** 2)):.2f} μGal") .. rst-class:: sphx-glr-script-out .. code-block:: none Gravity residuals: Mean: -0.00 μGal Std: 10540.46 μGal RMS: 10540.46 μGal .. GENERATED FROM PYTHON SOURCE LINES 437-450 Visualize the Model ------------------- Create 2D and 3D visualizations of the geological model **Model Visualization**: Before analyzing geophysical results, it's important to visualize the geological model itself to understand: * Geometry of geological units * Position of observation points relative to the modeled structures * Topographic influences * Overall model plausibility .. GENERATED FROM PYTHON SOURCE LINES 450-466 .. code-block:: Python import gempy_viewer as gpv import pyvista as pv # Configure PyVista to use smaller fonts for axes pv.global_theme.font.size = 10 pv.global_theme.font.label_size = 10 pv.global_theme.font.title_size = 12 gpv.plot_2d(simple_geo_model) plot_model_and_gravity_sensors(simple_geo_model, xy_ravel, grav) .. image-sg:: /examples_basic/images/sphx_glr_03_gravity_forward_model_001.png :alt: Cell Number: mid Direction: y :srcset: /examples_basic/images/sphx_glr_03_gravity_forward_model_001.png :class: sphx-glr-single-img .. image-sg:: /examples_basic/images/sphx_glr_03_gravity_forward_model_002.png :alt: 03 gravity forward model :srcset: /examples_basic/images/sphx_glr_03_gravity_forward_model_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/leguark/PycharmProjects/gempy_viewer/gempy_viewer/modules/plot_3d/drawer_surfaces_3d.py:38: PyVistaDeprecationWarning: ../../../gempy_viewer/gempy_viewer/modules/plot_3d/drawer_surfaces_3d.py:38: Argument 'color' must be passed as a keyword argument to function 'BasePlotter.add_mesh'. From version 0.50, passing this as a positional argument will result in a TypeError. gempy_vista.surface_actors[element.name] = gempy_vista.p.add_mesh( .. GENERATED FROM PYTHON SOURCE LINES 467-479 Visualize Gravity Results -------------------------- Plot forward gravity and comparison with observations **Geophysical Data Visualization**: Spatial plots help identify: * Anomaly patterns and their relationship to geology * Data coverage and survey design * Edge effects and boundary artifacts * Correlation between observed and modeled signatures .. GENERATED FROM PYTHON SOURCE LINES 479-491 .. code-block:: Python from mineye.GeoModel.plotting.probabilistic_analysis import plot_fw_geophysics, plot_comparison plot_fw_geophysics( fw_values=grav, observed_data=gravity_data, xy_ravel=xy_ravel, label=r'Forward Model Gravity ($\mu$gal)', title='Forward Model Gravity Results' ) plot_comparison(observed_gravity, grav, xy_ravel, unit_label=r'$\mu$Gal') .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_basic/images/sphx_glr_03_gravity_forward_model_003.png :alt: Forward Model Gravity Results :srcset: /examples_basic/images/sphx_glr_03_gravity_forward_model_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples_basic/images/sphx_glr_03_gravity_forward_model_004.png :alt: Observed (align_to_reference), Forward Model (align_to_reference), Residuals (Observed - Forward Model) (align_to_reference), Observed vs Forward Model Correlation :srcset: /examples_basic/images/sphx_glr_03_gravity_forward_model_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Plotting 260 actual measurement locations === Observed vs Predicted Comparison === Computing align_to_reference alignment parameters... Alignment params: {'method': 'align_to_reference', 'reference_mean': 25587.184615384616, 'reference_std': 8369.297471927344, 'baseline_forward_mean': -25587.18461538464, 'baseline_forward_std': 8353.187163491928} .. GENERATED FROM PYTHON SOURCE LINES 492-507 Visualize Forward Model Results -------------------------------- **Detailed Gravity Map**: This map shows the predicted gravity field from our geological model. Key features to observe: * **Positive anomalies**: Occur over the dense plutonite intrusion * **Negative anomalies**: May occur at the edges or over less dense rocks * **Spatial trends**: Reflect the geometry and extent of density contrasts * **Amplitude**: Related to the magnitude of density contrast and body volume The colormap (viridis_r) is reversed so higher values (positive anomalies) appear warmer, following common geophysical convention. .. GENERATED FROM PYTHON SOURCE LINES 507-532 .. code-block:: Python grav_forward = grav print(f"✓ Forward model computed") print(f" Gravity range: [{grav_forward.min():.2f}, {grav_forward.max():.2f}] µGal") fig, ax = plt.subplots(figsize=(10, 8)) scatter = ax.scatter( xy_ravel[:, 0], xy_ravel[:, 1], c=grav_forward, s=50, cmap='viridis_r', alpha=0.8, edgecolors='black', linewidth=0.5 ) cbar = plt.colorbar(scatter, ax=ax) cbar.set_label(r'Forward Model Gravity (µGal)', fontsize=12) ax.set_title('Forward Gravity Model', fontsize=14, fontweight='bold') ax.set_xlabel('X (m)') ax.set_ylabel('Y (m)') ax.grid(True, alpha=0.3) plt.tight_layout() plt.show() .. image-sg:: /examples_basic/images/sphx_glr_03_gravity_forward_model_005.png :alt: Forward Gravity Model :srcset: /examples_basic/images/sphx_glr_03_gravity_forward_model_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none ✓ Forward model computed Gravity range: [-50607.12, 29982.25] µGal .. GENERATED FROM PYTHON SOURCE LINES 533-561 Compare with Observed Gravity ------------------------------ **Three-Panel Comparison**: This comprehensive view shows: 1. **Observed Gravity** (Left): - Actual field measurements - Represents the "truth" we're trying to match - Contains geological signal + noise + regional trends 2. **Forward Model** (Center): - Predicted gravity from our geological interpretation - Shows what the model "thinks" the gravity should be - Ideally should match observed patterns 3. **Residuals** (Right): - Difference between observed and modeled (Obs - Model) - **Red**: Model under-predicts (observed > modeled) - **Blue**: Model over-predicts (observed < modeled) - **White/near-zero**: Good fit **What to look for**: * Systematic patterns in residuals suggest model improvements needed * Random residuals indicate we're at the noise level * Large residuals in specific areas point to local geological features not captured .. GENERATED FROM PYTHON SOURCE LINES 561-605 .. code-block:: Python # Convert observed from mGal to µGal observed_ugal = observed_gravity_ugal # Create comparison plot fig, axes = plt.subplots(1, 3, figsize=(18, 5)) # Observed sc1 = axes[0].scatter( xy_ravel[:, 0], xy_ravel[:, 1], c=observed_ugal, s=50, cmap='viridis_r', edgecolors='black', linewidth=0.5 ) axes[0].set_title('Observed Gravity') axes[0].set_xlabel('X (m)') axes[0].set_ylabel('Y (m)') plt.colorbar(sc1, ax=axes[0], label='µGal') # Forward model sc2 = axes[1].scatter( xy_ravel[:, 0], xy_ravel[:, 1], c=grav_forward, s=50, cmap='viridis_r', edgecolors='black', linewidth=0.5 ) axes[1].set_title('Forward Model') axes[1].set_xlabel('X (m)') plt.colorbar(sc2, ax=axes[1], label='µGal') # Residuals # Use a symmetric diverging colormap centered at zero res_max = np.max(np.abs(residuals)) sc3 = axes[2].scatter( xy_ravel[:, 0], xy_ravel[:, 1], c=residuals, s=50, cmap='RdBu_r', vmin=-res_max, vmax=res_max, edgecolors='black', linewidth=0.5 ) axes[2].set_title('Residuals (Obs - Model)') axes[2].set_xlabel('X (m)') plt.colorbar(sc3, ax=axes[2], label='µGal') plt.tight_layout() plt.show() .. image-sg:: /examples_basic/images/sphx_glr_03_gravity_forward_model_006.png :alt: Observed Gravity, Forward Model, Residuals (Obs - Model) :srcset: /examples_basic/images/sphx_glr_03_gravity_forward_model_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 606-630 Correlation Analysis -------------------- **Cross-Plot Analysis**: The scatter plot of observed vs. modeled gravity is a powerful diagnostic tool: * **1:1 line** (red dashed): Perfect agreement would have all points on this line * **Correlation coefficient (R)**: Measures linear relationship strength (-1 to +1) - R > 0.9: Excellent correlation - R = 0.7-0.9: Good correlation - R < 0.7: Poor correlation * **R²**: Fraction of variance explained by the model (0 to 1) * **RMSE**: Average prediction error in µGal **Interpreting the Plot**: * Points above the 1:1 line: Model under-predicts * Points below the 1:1 line: Model over-predicts * Scatter around the line: Combination of model error and measurement noise * Systematic deviations: Indicate model bias or missing physics This analysis helps determine if the geological model is a reasonable representation of the subsurface structure. .. GENERATED FROM PYTHON SOURCE LINES 630-663 .. code-block:: Python fig, ax = plt.subplots(figsize=(8, 8)) ax.scatter(observed_ugal, grav_forward, alpha=0.6, s=60, edgecolors='black', linewidth=0.5) # 1:1 line lims = [ min(ax.get_xlim()[0], ax.get_ylim()[0]), max(ax.get_xlim()[1], ax.get_ylim()[1]) ] ax.plot(lims, lims, 'r--', alpha=0.75, linewidth=2, label='1:1 line') # Calculate statistics correlation = np.corrcoef(observed_ugal, grav_forward)[0, 1] rmse = np.sqrt(np.mean(residuals ** 2)) ax.set_xlabel('Observed (µGal)', fontsize=12) ax.set_ylabel('Forward Model (µGal)', fontsize=12) ax.set_title('Observed vs Modeled Gravity', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) ax.legend() # Add statistics text box textstr = f'R = {correlation:.3f}\nR² = {correlation ** 2:.3f}\nRMSE = {rmse:.2f} µGal' ax.text(0.05, 0.95, textstr, transform=ax.transAxes, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8), fontsize=11) plt.tight_layout() plt.show() .. image-sg:: /examples_basic/images/sphx_glr_03_gravity_forward_model_007.png :alt: Observed vs Modeled Gravity :srcset: /examples_basic/images/sphx_glr_03_gravity_forward_model_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 664-666 Print Summary Statistics ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 666-687 .. code-block:: Python print("\n" + "=" * 50) print("GRAVITY COMPARISON STATISTICS") print("=" * 50) print(f"Number of measurements: {len(observed_ugal)}") print(f"\nObserved gravity:") print(f" Mean: {observed_ugal.mean():.2f} µGal") print(f" Std: {observed_ugal.std():.2f} µGal") print(f" Range: [{observed_ugal.min():.2f}, {observed_ugal.max():.2f}] µGal") print(f"\nForward model:") print(f" Mean: {grav_forward.mean():.2f} µGal") print(f" Std: {grav_forward.std():.2f} µGal") print(f" Range: [{grav_forward.min():.2f}, {grav_forward.max():.2f}] µGal") print(f"\nResiduals:") print(f" Mean: {residuals.mean():.2f} µGal") print(f" Std: {residuals.std():.2f} µGal") print(f" RMSE: {rmse:.2f} µGal") print(f" MAE: {np.abs(residuals).mean():.2f} µGal") print(f"\nCorrelation: {correlation:.4f} (R² = {correlation ** 2:.4f})") print("=" * 50) .. rst-class:: sphx-glr-script-out .. code-block:: none ================================================== GRAVITY COMPARISON STATISTICS ================================================== Number of measurements: 260 Observed gravity: Mean: 25587.18 µGal Std: 8369.30 µGal Range: [9508.00, 47323.00] µGal Forward model: Mean: 25587.18 µGal Std: 8353.19 µGal Range: [-50607.12, 29982.25] µGal Residuals: Mean: -0.00 µGal Std: 10540.46 µGal RMSE: 10540.46 µGal MAE: 7355.18 µGal Correlation: 0.2054 (R² = 0.0422) ================================================== .. GENERATED FROM PYTHON SOURCE LINES 688-725 Summary and Interpretation --------------------------- **Model Limitations**: Forward models are simplifications of reality and have inherent limitations: * **Density assumptions**: We used uniform densities, but real rocks are heterogeneous * **Geological simplification**: Simple geometry may not capture all structural complexity * **Data coverage**: Limited measurement points constrain resolution * **Regional trends**: Long-wavelength features may not be fully captured * **3D effects**: Some anomalies may be caused by out-of-plane structures **Next Steps and Applications**: This forward modeling foundation enables: * **Probabilistic Inversion** (Example 04): Use Bayesian methods to infer densities and geological parameters from the gravity data * **Uncertainty Quantification**: Propagate uncertainties in geology to geophysics * **Joint Inversion**: Combine gravity with magnetics, seismic, or other data * **Model Refinement**: Use residuals to iteratively improve geological interpretation * **Resource Estimation**: Better constrain ore body geometry and properties **Further Reading**: * Blakely, R. J. (1996). Potential Theory in Gravity and Magnetic Applications. Cambridge University Press. * Li, Y., & Oldenburg, D. W. (1998). 3-D inversion of gravity data. Geophysics, 63(1), 109-119. * de la Varga et al. (2019). GemPy 1.0: open-source stochastic geological modeling and inversion. Geoscientific Model Development, 12(1), 1-32. .. seealso:: * :doc:`../examples_probabilistic/04_gravity_inversion` - Bayesian gravity inversion * :doc:`../examples_probabilistic/05_magnetics_inversion` - Joint geophysical inversion sphinx_gallery_thumbnail_number = 2 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.562 seconds) .. _sphx_glr_download_examples_basic_03_gravity_forward_model.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03_gravity_forward_model.ipynb <03_gravity_forward_model.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03_gravity_forward_model.py <03_gravity_forward_model.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03_gravity_forward_model.zip <03_gravity_forward_model.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_