{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calibrated and not calibrated ERT with Boxford dataset\n",
    "This study has been demonstrated in the EMagPy paper [McLachlan et al. (2021)](https://doi.org/10.1016/j.cageo.2020.104561).\n",
    "\n",
    "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)).\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import sys\n",
    "sys.path.append('../src') # add path where emagpy is\n",
    "from emagpy import Problem\n",
    "\n",
    "datadir = '../src/examples/'\n",
    "\n",
    "letters = ['a','b','c','d','e','f','g','h','i','j']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fnameEC = datadir + 'boxford-calib/eri_ec.csv'\n",
    "fnameECa = datadir + 'boxford-calib/eca_calibration.csv'\n",
    "\n",
    "# non calibrated\n",
    "k1 = Problem()\n",
    "k1.createSurvey(fnameECa)\n",
    "k1.show()\n",
    "k1.setInit(depths0=np.arange(0.05, 3, 0.25))\n",
    "k1.invert(forwardModel='FSlin', alpha=0.01, method='L-BFGS-B', njobs=-1)\n",
    "\n",
    "# ERT calibrated\n",
    "k2 = Problem()\n",
    "k2.createSurvey(fnameECa)\n",
    "k2.calibrate(fnameECa, fnameEC, forwardModel='FSlin') # plot calibration\n",
    "k2.calibrate(fnameECa, fnameEC, forwardModel='FSlin', apply=True) # apply the calibration\n",
    "k2.show()\n",
    "k2.setInit(depths0=np.arange(0.05, 3, 0.25))\n",
    "k2.invert(forwardModel='FSlin', alpha=0.001, method='L-BFGS-B', njobs=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# figure\n",
    "fig, axs = plt.subplots(3, 1, sharey=True, sharex=True, figsize=(10,6))\n",
    "peatdepths = pd.read_csv(datadir + 'boxford-calib/peat-depth.dat', sep=\"\\t\")\n",
    "ax = axs[0]\n",
    "k1.showResults(ax=ax, vmin=0, vmax=35)\n",
    "ax.plot(peatdepths['distance (m)'], -peatdepths['depth (m)'], 'r--')\n",
    "ax.set_title('(a) Uncalibrated')\n",
    "ax.set_xlabel('')\n",
    "ax = axs[1]\n",
    "k2.showResults(ax=ax, vmin=0, vmax=35)\n",
    "ax.plot(peatdepths['distance (m)'], -peatdepths['depth (m)'], 'r--')\n",
    "ax.set_title('(b) ERT calibrated')\n",
    "ax.set_xlabel('')\n",
    "ax = axs[2]\n",
    "kres = Problem()\n",
    "kres.importModel(datadir + 'boxford-calib/ert-section.csv')\n",
    "kres.models[0] = kres.models[0][::4,:] # because taken every 0.25 m\n",
    "kres.depths[0] = kres.depths[0][::4,:]\n",
    "kres.showResults(ax=ax, vmin=0, vmax=35, maxDepth=1.5)\n",
    "ax.set_title('(c) ERT transect')\n",
    "ax.plot(peatdepths['distance (m)'], -peatdepths['depth (m)'], 'r--')\n",
    "ax.set_xlim([0, 40])\n",
    "ax.set_ylim([-3, 0])\n",
    "ax.set_xlabel('Distance [m]')\n",
    "ax.set_aspect('auto')\n",
    "fig.tight_layout()\n",
    "fig.savefig(outputdir + 'fig10calibrated-inv.jpg', dpi=500)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}