.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/Ex_gravity_magnetics_inversion.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_auto_examples_Ex_gravity_magnetics_inversion.py: Ex. Gravity and Magnetics Processing and Inversion ================================================== This example uses the promoted, Qt-free potential-field API to inspect field datasets, separate regional and residual anomalies, evaluate analytic body responses, and run compact 3D SimPEG inversions. .. GENERATED FROM PYTHON SOURCE LINES 9-34 .. code-block:: Python from pathlib import Path import sys current_dir = Path(__file__).resolve().parent if "__file__" in globals() else Path.cwd() repo_root = current_dir.parent if str(repo_root) not in sys.path: sys.path.insert(0, str(repo_root)) import matplotlib.pyplot as plt import numpy as np import pandas as pd from PyHydroGeophysX.workflows.gravmag import ( forward_bodies, grid_data, invert_gravmag, regional_residual, ) data_dir = current_dir / "data" / "Gravity_Magnetics" output_dir = current_dir / "results" / "gravity_magnetics_inversion" output_dir.mkdir(parents=True, exist_ok=True) .. GENERATED FROM PYTHON SOURCE LINES 36-38 Load and QC the Field Data -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 40-71 .. code-block:: Python gravity = pd.read_csv(data_dir / "bushveld_gravity_disturbance.csv") magnetic = pd.read_csv(data_dir / "britain_aeromagnetic_anomaly.csv") gx = gravity.iloc[:, 0].to_numpy(float) gy = gravity.iloc[:, 1].to_numpy(float) gobs = gravity.iloc[:, 2].to_numpy(float) mx = magnetic.iloc[:, 0].to_numpy(float) my = magnetic.iloc[:, 1].to_numpy(float) mobs = magnetic.iloc[:, 2].to_numpy(float) gregional, gresidual = regional_residual(gx, gy, gobs, degree=1) mregional, mresidual = regional_residual(mx, my, mobs, degree=1) ggrid = grid_data(gx, gy, gresidual, nx=100, ny=100) mgrid = grid_data(mx, my, mresidual, nx=100, ny=100) fig, axes = plt.subplots(2, 2, figsize=(12, 8), constrained_layout=True) for ax, x, y, values, title, unit in ( (axes[0, 0], gx, gy, gobs, "Bushveld observed gravity", "mGal"), (axes[0, 1], gx, gy, gresidual, "Bushveld residual gravity", "mGal"), (axes[1, 0], mx, my, mobs, "Britain observed magnetics", "nT"), (axes[1, 1], mx, my, mresidual, "Britain residual magnetics", "nT"), ): artist = ax.scatter(x, y, c=values, s=5, cmap="RdBu_r") fig.colorbar(artist, ax=ax, label=unit) ax.set_title(title) ax.set_xlabel("Easting (m)") ax.set_ylabel("Northing (m)") ax.set_aspect("equal", adjustable="box") fig.savefig(output_dir / "Ex_gravity_magnetics_inversion_fig_01.png", dpi=240) plt.show() .. GENERATED FROM PYTHON SOURCE LINES 72-78 The paired maps distinguish broad regional trends from shorter-wavelength anomalies that are appropriate inputs for compact-source interpretation. .. image:: /auto_examples/images/Ex_gravity_magnetics_inversion_fig_01.png :align: center :width: 800px .. GENERATED FROM PYTHON SOURCE LINES 80-82 Analytic Forward Responses -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 84-110 .. code-block:: Python xline = np.linspace(-200.0, 200.0, 301) yline = np.zeros_like(xline) gravity_body = [{ "type": "sphere", "x0": 0.0, "y0": 0.0, "z0": 80.0, "radius": 35.0, "density_contrast": 450.0, }] magnetic_body = [{ "type": "sphere", "x0": 30.0, "y0": 0.0, "z0": 90.0, "radius": 30.0, "susceptibility": 0.03, }] field = {"strength": 50000.0, "inclination": 65.0, "declination": 0.0} gforward = forward_bodies(xline, yline, "gravity", gravity_body) mforward = forward_bodies(xline, yline, "magnetics", magnetic_body, field=field) fig, axes = plt.subplots(1, 2, figsize=(11, 4), constrained_layout=True) axes[0].plot(xline, gforward, color="tab:blue") axes[0].set(xlabel="Distance (m)", ylabel="Gravity anomaly (mGal)", title="Buried sphere gravity response") axes[1].plot(xline, mforward, color="tab:red") axes[1].set(xlabel="Distance (m)", ylabel="TMI anomaly (nT)", title="Induced magnetic response") for ax in axes: ax.grid(alpha=0.25) fig.savefig(output_dir / "Ex_gravity_magnetics_inversion_fig_02.png", dpi=240) plt.show() .. GENERATED FROM PYTHON SOURCE LINES 111-114 .. image:: /auto_examples/images/Ex_gravity_magnetics_inversion_fig_02.png :align: center :width: 800px .. GENERATED FROM PYTHON SOURCE LINES 116-121 Compact 3D Field-Data Inversions -------------------------------- Small meshes keep this tutorial practical. The recovered density and susceptibility contrasts should be interpreted qualitatively. .. GENERATED FROM PYTHON SOURCE LINES 123-159 .. code-block:: Python gresult = invert_gravmag( gx, gy, gobs, "gravity", detrend=1, max_stations=250, n_xy=8, n_z=4, max_iterations=4, ) mresult = invert_gravmag( mx, my, mobs, "magnetics", detrend=1, max_stations=220, n_xy=7, n_z=4, max_iterations=4, field={"strength_nT": 48800.0, "inclination": 66.0, "declination": -2.0}, ) fig, axes = plt.subplots(2, 2, figsize=(12, 8), constrained_layout=True) for row, result, title in ( (0, gresult, "Gravity density contrast"), (1, mresult, "Magnetic susceptibility contrast"), ): ex, ey, _ = result["edges"] model = np.asarray(result["model3d"]) image = axes[row, 0].pcolormesh(ex, ey, model[:, :, -1].T, shading="auto", cmap=result["cmap"]) fig.colorbar(image, ax=axes[row, 0], label=result["label"]) axes[row, 0].set_title(f"{title}: shallowest layer") axes[row, 0].set_xlabel("Easting (m)") axes[row, 0].set_ylabel("Northing (m)") history = np.asarray(result.get("convergence", []), dtype=float) axes[row, 1].plot(np.arange(1, history.size + 1), history, "o-") axes[row, 1].axhline(1.0, color="k", ls="--", lw=1) axes[row, 1].set_title(f"Convergence (final chi²={result['chi2']:.2f})") axes[row, 1].set_xlabel("Iteration") axes[row, 1].set_ylabel("Normalized misfit") axes[row, 1].grid(alpha=0.25) fig.savefig(output_dir / "Ex_gravity_magnetics_inversion_fig_03.png", dpi=240) plt.show() print(f"Gravity inversion: {gresult['n_cells']} cells, chi2={gresult['chi2']:.3f}") print(f"Magnetic inversion: {mresult['n_cells']} cells, chi2={mresult['chi2']:.3f}") .. GENERATED FROM PYTHON SOURCE LINES 160-163 .. image:: /auto_examples/images/Ex_gravity_magnetics_inversion_fig_03.png :align: center :width: 800px .. _sphx_glr_download_auto_examples_Ex_gravity_magnetics_inversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Ex_gravity_magnetics_inversion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Ex_gravity_magnetics_inversion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Ex_gravity_magnetics_inversion.zip `