HLEM Ground based Maxmin 100m cable with topography

[ajelenic]

Question

I saw that there is a method to import Apex Maxmin in-phase %, and quadrature % data, but there are no examples that I could find.
Here is my attempt, which gives the error "inversion failed: 'FDEM' object has no attribute 'invert'"

import numpy as np
import pandas as pd
import pygimli as pg
from pygimli.physics.em import fdem
import matplotlib.pyplot as plt

# -------------------------------
# 1. Load MaxMin CSV
# -------------------------------
csv_path = r"input\L5850.csv"
df = pd.read_csv(csv_path)

# Frequency columns
freq_cols_in = ["IP_880", "IP_3520", "IP_14080"]
freq_cols_q  = ["Q_880",  "Q_3520",  "Q_14080"]

# Frequencies for the forward operator (must match number of columns)
frequencies = pg.Vector([880.0, 3520.0, 14080.0])  # Hz

# Station coordinates
x = df["X"].values

# Data arrays
inphase = df[freq_cols_in].to_numpy(dtype=float)
quad    = df[freq_cols_q].to_numpy(dtype=float)

# Mask invalid stations
mask_valid = ~np.isnan(inphase + quad).any(axis=1)
print(f"Valid stations: {mask_valid.sum()} / {len(mask_valid)}")

xv = x[mask_valid]
inphase_v = inphase[mask_valid]
quad_v    = quad[mask_valid]

n_stations = len(xv)

# ---------------------------------------------------------------------
# 2. Run inversion per station
# ---------------------------------------------------------------------
nlay = 3
coilspacing = 100.0  # meters
z = 1.25            # measurement height

inversion_results = np.full((n_stations, nlay), np.nan)

for i in range(n_stations):
    try:
        # Create FDEM 1D forward operator for this station
        fop = fdem.FDEM(
            nlay,
            frequencies,
            coilspacing,
            z,
          
        )

        # Prepare complex data (in-phase + i*quadrature)
        d_complex = inphase_v[i, :] + 1j * quad_v[i, :]

        # For pyGIMLi inversion, concatenate real + imag parts
        inv_data = pg.Vector(np.hstack((d_complex.real, d_complex.imag)))

        # Initial model (S/m)
        start_model = pg.Vector([0.01] * nlay)

        # Run 1D inversion
        model = fop.invert(inv_data, start_model)
        inversion_results[i, :] = model

    except Exception as e:
        print(f"Station {i} inversion failed: {e}")

# ---------------------------------------------------------------------
# 3. Plot section
# ---------------------------------------------------------------------
layer_thicknesses = [10, 20, 50]
depth_mid = np.cumsum(layer_thicknesses) - np.array(layer_thicknesses)/2

plt.figure(figsize=(12, 6))
plt.pcolormesh(xv, depth_mid, inversion_results.T,
               cmap="plasma", shading="auto")
plt.gca().invert_yaxis()
plt.xlabel("X [m]")
plt.ylabel("Depth [m]")
plt.colorbar(label="Conductivity [S/m]")
plt.title("1D FDEM Inversion (pyGIMLi FDEM1dModelling, Line L5850)")
plt.show()

[halbmy]

There is no such method as fop.invert(). You need to setup an inversion instance, feed the forward operator, a data and error vector and run it.

[ajelenic]

I have updated my script, can you look at the fop.
TypeError: unsupported operand type(s) for *: 'float' and 'NoneType'
I think this is related to my misunderstanding that r is related to coil seperation.

frequencies = [880, 3520, 14080]
coil_spacing =99.9
nlay=5
sensor_height =1.0
coil_spacing =99.9

fdem_data.fop = HEMmodelling(
frequencies=frequencies, # must be first
nlay=nlay, # scalar
height=sensor_height, # scalar
coil_spacing=coil_spacing, # scalar
r=coil_spacing

)