#!/usr/bin/env python
# coding: utf-8

# ## Calibrated and not calibrated ERT with Boxford dataset
# This study has been demonstrated in the EMagPy paper [McLachlan et al. (2021)](https://doi.org/10.1016/j.cageo.2020.104561).
# 
# 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](https://doi.org/10.3997/1873-0604.2010037)).
# 
# 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.

# In[ ]:


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']


# In[ ]:


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)


# In[ ]:


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