.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/Ex_EM_line_section.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_EM_line_section.py: Ex. Airborne EM Line Inversion ============================== This example loads a VTEM line and its survey geometry, inspects measured decays, calibrates the response to a documented reference resistivity, and creates a stitched resistivity section from independent 1D inversions. .. GENERATED FROM PYTHON SOURCE LINES 9-37 .. 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 from matplotlib.colors import LogNorm import numpy as np import pandas as pd from PyHydroGeophysX.workflows.em1d import ( calibrate_to_reference, invert_line, load_line_geometry, load_sounding, ) data_dir = current_dir / "data" / "EM" / "EastRiver_VTEM" output_dir = current_dir / "results" / "em_line_section" output_dir.mkdir(parents=True, exist_ok=True) data_path = data_dir / "eastriver_vtem_line22030.csv" geometry_path = data_dir / "eastriver_vtem_line22030_geometry.csv" .. GENERATED FROM PYTHON SOURCE LINES 39-41 Inspect the VTEM Soundings -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 43-60 .. code-block:: Python geometry = load_line_geometry(str(geometry_path)) fig, ax = plt.subplots(figsize=(8, 5), constrained_layout=True) for sounding in (0, 5, 10, 15, 21): data = load_sounding(str(data_path), "TDEM", sounding=sounding) ax.loglog(data["times"], np.abs(data["response"]), marker="o", ms=3, label=f"Sounding {sounding + 1}") ax.set_xlabel("Time (s)") ax.set_ylabel("|Response|") ax.set_title("East River VTEM decay curves") ax.grid(True, which="both", alpha=0.25) ax.legend(ncols=2) fig.savefig(output_dir / "Ex_EM_line_section_fig_01.png", dpi=240) plt.show() print(f"Loaded {load_sounding(str(data_path), 'TDEM')['n_soundings']} soundings") print(f"Line distance: {geometry['positions'].min():.1f}–{geometry['positions'].max():.1f} m") .. GENERATED FROM PYTHON SOURCE LINES 61-64 .. image:: /auto_examples/images/Ex_EM_line_section_fig_01.png :align: center :width: 700px .. GENERATED FROM PYTHON SOURCE LINES 66-71 Calibrate and Invert the Line ----------------------------- The documented 320 ohm-m reference constrains the otherwise ambiguous amplitude scale. Twelve soundings are used here to keep the tutorial compact. .. GENERATED FROM PYTHON SOURCE LINES 73-123 .. code-block:: Python forward_geometry = { "source_radius": 13.0, "height": 82.0, "orientation": "z", "waveform": "step_off", } inversion_parameters = { "n_layers": 9, "min_thickness": 8.0, "max_thickness": 55.0, "starting_resistivity": 320.0, "smoothness": 0.5, "rel_error": 0.08, "max_iterations": 6, } data_scale = calibrate_to_reference( str(data_path), "TDEM", forward_geometry, inversion_parameters, ref_resistivity=320.0, max_probe=4, ) inversion_parameters["data_scale"] = data_scale result = invert_line( str(data_path), "TDEM", forward_geometry, inversion_parameters, positions=geometry["positions"], heights=geometry["heights"], max_soundings=12, doi_blank=True, out_dir=output_dir, ) ex, _, ez = result["edges"] section = np.asarray(result["model3d"])[:, 0, :].T finite = section[np.isfinite(section) & (section > 0)] vmin, vmax = np.nanpercentile(finite, [5, 95]) fig, axes = plt.subplots(2, 1, figsize=(11, 7), constrained_layout=True, gridspec_kw={"height_ratios": [3, 1]}) image = axes[0].pcolormesh(ex, ez, section, shading="auto", cmap="turbo", norm=LogNorm(vmin=max(vmin, 1.0), vmax=vmax)) fig.colorbar(image, ax=axes[0], label="Resistivity (ohm m)") axes[0].set_ylabel("Elevation relative to surface (m)") axes[0].set_title("Independent 1D VTEM inversions stitched along line") chi2 = np.asarray(result["chi2_list"], dtype=float) axes[1].bar(result["positions"], chi2, width=25.0, color="tab:blue") axes[1].axhline(1.0, color="k", ls="--", lw=1) axes[1].set_xlabel("Distance along line (m)") axes[1].set_ylabel("chi²") axes[1].grid(axis="y", alpha=0.25) fig.savefig(output_dir / "Ex_EM_line_section_fig_02.png", dpi=240) plt.show() print(f"Calibration scale: {data_scale:.4g}") print(f"Mean line chi2: {result['chi2']:.3f}") .. GENERATED FROM PYTHON SOURCE LINES 124-127 .. image:: /auto_examples/images/Ex_EM_line_section_fig_02.png :align: center :width: 800px .. GENERATED FROM PYTHON SOURCE LINES 129-144 Compare with the Published Reference Section -------------------------------------------- The comparison is restricted to the shallowest 100 m, which is what this data set can support. The recorded gates stop at 391 microseconds, giving a diffusion depth near 450 m in a 320 ohm-m half-space, and sensitivity has already fallen off well above that. The published USGS section extends to 1000 m; below roughly 100 m the inversion here drifts progressively resistive (a factor of four by 200 m) because the late-time data that would hold it down were never recorded. Comparing the full depth range would say more about that missing sensitivity than about either inversion. The amplitude is a second reason to read this as a check rather than a validation: ``calibrate_to_reference`` sets the data scale by assuming a 320 ohm-m half-space, so it is an assumption, not an instrument calibration. .. GENERATED FROM PYTHON SOURCE LINES 146-218 .. code-block:: Python COMPARISON_DEPTH_M = 100.0 reference = pd.read_csv(data_dir / "eastriver_vtem_line22030_usgs_resistivity.csv") reference_depth = reference.iloc[:, 0].to_numpy(float) reference_model = reference.iloc[:, 1:13].to_numpy(float) positions = np.asarray(result["positions"], dtype=float) depth_center = -0.5 * (ez[:-1] + ez[1:]) # Both grids ordered by increasing depth, then cut at the comparison depth. order = np.argsort(depth_center) depth_center, section_sorted = depth_center[order], section[order] ours_band = depth_center <= COMPARISON_DEPTH_M ref_band = reference_depth <= COMPARISON_DEPTH_M # One colour scale for both panels, taken from both panels. Scaling it to this # inversion alone (as a percentile of `section`) pushes the reference into the # bottom of the colour map and hides its structure. paired = np.concatenate([ section_sorted[ours_band].ravel(), reference_model[ref_band].ravel()]) paired = paired[np.isfinite(paired) & (paired > 0)] norm = LogNorm(*np.nanpercentile(paired, [2, 98])) fig, axes = plt.subplots(1, 2, figsize=(12, 5), constrained_layout=True, sharey=True) axes[0].pcolormesh(positions, depth_center[ours_band], section_sorted[ours_band], shading="nearest", cmap="turbo", norm=norm) axes[0].set_title("PyHydroGeophysX") axes[0].set_xlabel("Distance (m)") axes[0].set_ylabel("Depth (m)") ref_image = axes[1].pcolormesh(positions, reference_depth[ref_band], reference_model[ref_band], shading="nearest", cmap="turbo", norm=norm) axes[1].set_title("Provided USGS reference") axes[1].set_xlabel("Distance (m)") axes[0].set_ylim(COMPARISON_DEPTH_M, 0) fig.colorbar(ref_image, ax=axes, label="Resistivity (ohm m)") fig.suptitle(f"Upper {COMPARISON_DEPTH_M:.0f} m, the depth range these gates constrain") fig.savefig(output_dir / "Ex_EM_line_section_fig_03.png", dpi=240) plt.show() # Score the overlap instead of leaving the reader to judge two colour maps. on_reference = np.full((ref_band.sum(), section_sorted.shape[1]), np.nan) for column in range(section_sorted.shape[1]): trace = section_sorted[:, column] usable = np.isfinite(trace) & (trace > 0) if usable.sum() >= 2: on_reference[:, column] = np.interp( reference_depth[ref_band], depth_center[usable], trace[usable], left=np.nan, right=np.nan) comparable = np.isfinite(on_reference) & np.isfinite(reference_model[ref_band]) ours_log = np.log10(on_reference[comparable]) theirs_log = np.log10(reference_model[ref_band][comparable]) print(f"Upper {COMPARISON_DEPTH_M:.0f} m, {comparable.sum()} comparable cells") print(f" log10 RMSE vs USGS: {np.sqrt(np.mean((ours_log - theirs_log) ** 2)):.3f}" f" (0.30 is a factor of two)") print(f" median resistivity: {10 ** np.median(ours_log):.0f} ohm m here, " f"{10 ** np.median(theirs_log):.0f} ohm m in the reference") # Per band, so the depth at which agreement decays is visible as a number. print(" by depth band:") for lower, upper in ((0.0, 40.0), (40.0, 70.0), (70.0, COMPARISON_DEPTH_M)): rows = (reference_depth[ref_band] >= lower) & (reference_depth[ref_band] < upper) cells = comparable[rows] if cells.sum() < 3: continue band_ours = np.log10(on_reference[rows][cells]) band_theirs = np.log10(reference_model[ref_band][rows][cells]) print(f" {lower:>3.0f}-{upper:>3.0f} m: " f"log10 RMSE {np.sqrt(np.mean((band_ours - band_theirs) ** 2)):.3f}, " f"median {10 ** np.median(band_ours):>4.0f} vs " f"{10 ** np.median(band_theirs):.0f} ohm m") .. GENERATED FROM PYTHON SOURCE LINES 219-234 .. image:: /auto_examples/images/Ex_EM_line_section_fig_03.png :align: center :width: 800px Agreement is closest between about 40 and 70 m, where the two sections differ by a factor of 1.5 (log10 RMSE 0.167, medians 211 and 257 ohm-m). Over the whole upper 100 m it is a factor of two (log10 RMSE 0.274), which is what to expect from independent 1D inversions using different codes, layer grids and regularization. Both sections also place the same near-surface conductor around 1200 to 1400 m along the line. The deepest row drawn here is the layer centred at 91 m, and it is already running high: 523 ohm-m against 280 in the reference. That row is the start of the resistive drift described above, so read the bottom of this figure as the edge of the usable range rather than as a result. .. _sphx_glr_download_auto_examples_Ex_EM_line_section.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Ex_EM_line_section.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Ex_EM_line_section.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Ex_EM_line_section.zip `