Calibrated and not calibrated ERT with Boxford dataset

This study has been demonstrated in the EMagPy paper McLachlan et al. (2021).

To obtain quantitative EMI measurements, a calibration is needed. One way to perform this calibration is to use an inverted resistivity model from ERT and perform a forward EM model on it. We then match the generated ECa with the ECa measured on the ERT transect (see Lavoué et al. 2010).

The results are shown in the figure below. Smoothly inverted non-calibrated (a) and calibrated (b) EMI data with the corresponding ERT inversion (c). The red line shows the true depth of the peat intrusive penetration measurements.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sys
sys.path.append('../src') # add path where emagpy is
from emagpy import Problem

datadir = '../src/examples/'

letters = ['a','b','c','d','e','f','g','h','i','j']

fnameEC = datadir + 'boxford-calib/eri_ec.csv'
fnameECa = datadir + 'boxford-calib/eca_calibration.csv'

# non calibrated
k1 = Problem()
k1.createSurvey(fnameECa)
k1.show()
k1.setInit(depths0=np.arange(0.05, 3, 0.25))
k1.invert(forwardModel='FSlin', alpha=0.01, method='L-BFGS-B', njobs=-1)

# ERT calibrated
k2 = Problem()
k2.createSurvey(fnameECa)
k2.calibrate(fnameECa, fnameEC, forwardModel='FSlin') # plot calibration
k2.calibrate(fnameECa, fnameEC, forwardModel='FSlin', apply=True) # apply the calibration
k2.show()
k2.setInit(depths0=np.arange(0.05, 3, 0.25))
k2.invert(forwardModel='FSlin', alpha=0.001, method='L-BFGS-B', njobs=-1)

0%|          | 0/43 [00:00<?, ?it/s]

Survey 1/1

100%|██████████| 43/43 [00:16<00:00,  2.66it/s]

VCP1.48f10000h1: ECa(ERT) = 0.25 * ECa(EMI) +2.07 (R^2=0.38)
VCP2.82f10000h1: ECa(ERT) = 0.40 * ECa(EMI) +2.33 (R^2=0.48)
VCP4.49f10000h1: ECa(ERT) = 0.48 * ECa(EMI) +1.86 (R^2=0.53)
HCP1.48f10000h1: ECa(ERT) = 0.49 * ECa(EMI) +3.41 (R^2=0.36)
HCP2.82f10000h1: ECa(ERT) = 0.73 * ECa(EMI) +1.55 (R^2=0.54)
HCP4.49f10000h1: ECa(ERT) = 0.59 * ECa(EMI) +1.85 (R^2=0.36)

0%|          | 0/43 [00:00<?, ?it/s]

VCP1.48f10000h1: ECa(ERT) = 0.25 * ECa(EMI) +2.07 (R^2=0.38)
VCP2.82f10000h1: ECa(ERT) = 0.40 * ECa(EMI) +2.33 (R^2=0.48)
VCP4.49f10000h1: ECa(ERT) = 0.48 * ECa(EMI) +1.86 (R^2=0.53)
HCP1.48f10000h1: ECa(ERT) = 0.49 * ECa(EMI) +3.41 (R^2=0.36)
HCP2.82f10000h1: ECa(ERT) = 0.73 * ECa(EMI) +1.55 (R^2=0.54)
HCP4.49f10000h1: ECa(ERT) = 0.59 * ECa(EMI) +1.85 (R^2=0.36)
Correction is applied.
Survey 1/1

100%|██████████| 43/43 [00:45<00:00,  1.07s/it]
../_images/nb_paper-ert-calibration_0.png ../_images/nb_paper-ert-calibration_1.png ../_images/nb_paper-ert-calibration_2.png ../_images/nb_paper-ert-calibration_3.png 
# figure
fig, axs = plt.subplots(3, 1, sharey=True, sharex=True, figsize=(10,6))
peatdepths = pd.read_csv(datadir + 'boxford-calib/peat-depth.dat', sep="\t")
ax = axs[0]
k1.showResults(ax=ax, vmin=0, vmax=35)
ax.plot(peatdepths['distance (m)'], -peatdepths['depth (m)'], 'r--')
ax.set_title('(a) Uncalibrated')
ax.set_xlabel('')
ax = axs[1]
k2.showResults(ax=ax, vmin=0, vmax=35)
ax.plot(peatdepths['distance (m)'], -peatdepths['depth (m)'], 'r--')
ax.set_title('(b) ERT calibrated')
ax.set_xlabel('')
ax = axs[2]
kres = Problem()
kres.importModel(datadir + 'boxford-calib/ert-section.csv')
kres.models[0] = kres.models[0][::4,:] # because taken every 0.25 m
kres.depths[0] = kres.depths[0][::4,:]
kres.showResults(ax=ax, vmin=0, vmax=35, maxDepth=1.5)
ax.set_title('(c) ERT transect')
ax.plot(peatdepths['distance (m)'], -peatdepths['depth (m)'], 'r--')
ax.set_xlim([0, 40])
ax.set_ylim([-3, 0])
ax.set_xlabel('Distance [m]')
ax.set_aspect('auto')
fig.tight_layout()
fig.savefig(outputdir + 'fig10calibrated-inv.jpg', dpi=500)
../_images/nb_paper-ert-calibration_4.png

Download python script: nb_paper-ert-calibration.py

Download Jupyter notebook: nb_paper-ert-calibration.ipynb

View the notebook in the Jupyter nbviewer

Run this example interactively: binder