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.

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)

Load and QC the Field Data

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()

The paired maps distinguish broad regional trends from shorter-wavelength anomalies that are appropriate inputs for compact-source interpretation.

Analytic Forward Responses

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()

Compact 3D Field-Data Inversions

Small meshes keep this tutorial practical. The recovered density and susceptibility contrasts should be interpreted qualitatively.

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}")