{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Checkout www.pygimli.org for more examples"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Incorporating prior data into ERT inversion\n\nPrior data can often help to overcome ambiguity in the inversion process.\nHere we demonstrate the use of direct push (DP) data in an ERT inversion of\ndata collected at the surface.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nimport pygimli as pg\nfrom pygimli.physics import ert\nfrom pygimli.frameworks import PriorModelling, JointModelling\nfrom pygimli.viewer.mpl import draw1DColumn"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The prior data\n\nThis field data is from a site with layered sands and clays over a\nresistive bedrock. We load it from the example repository.\n\nAs a position of x=155m (center of the profile) we have a\nborehole/direct push with known in-situ-data. We load the three-column\nfile using numpy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x, z, r = pg.getExampleData(\"ert/bedrock.txt\").T\nfig, ax = plt.subplots()\nax.semilogx(r, z, \"*-\")\nax.set_xlabel(r\"$\\rho$ ($\\Omega$m)\")\nax.set_ylabel(\"depth (m)\")\nax.grid(True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We mainly see four layers: 1. a conductive (clayey) overburden of about\n17m thickness, 2. a medium resistivity interbedding of silt and sand,\nabout 7m thick 3. again clay with 8m thickness 4. the resistive bedrock\nwith a few hundred $\\Omega$ m\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The ERT data\nWe load the ERT data from the example repository and plot the pseudosection.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data = pg.getExampleData(\"ert/bedrock.dat\")\nprint(data)\nax, cb = ert.show(data)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The apparent resistivities show increasing values with larger spacings\nwith no observable noise. We first compute geometric factors and\nestimate an error model using rather low values for both error parts.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data.estimateError(relativeError=0.025, absoluteUError=100e-6)\n# data[\"err\"] = ert.estimateError(data, relativeError=0.025, absoluteUError=100e-6)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create an ERT manager and invert the data, already using a rather low value for\nthe vertical smoothness to account for the layered model.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mgr = ert.ERTManager(data, verbose=True)\nmgr.invert(paraDepth=70, quality=34, paraMaxCellSize=3000, zWeight=0.15, lam=30)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For reasons of comparability, we define a unique colormap and store all\noptions in a dictionary to be used in subsequent show commands.\n\nWe plot the result with these and plot the DP points onto the mesh.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kw = dict(cMin=20, cMax=500, logScale=True, cMap=\"Spectral_r\",\n          xlabel=\"x (m)\", ylabel=\"y (m)\")\nax, cb = mgr.showResult(**kw)\nzz = np.abs(z)\niz = np.argsort(z)\ndz = np.diff(zz[iz])\nthk = np.hstack([dz, dz[-1]])\nztop = -zz[iz[0]]-dz[0]/2\ncolkw = dict(x=x[0], val=r[iz], thk=thk, width=4, ztopo=ztop)\ndraw1DColumn(ax, **colkw, **kw)\nax.grid(True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We want to extract the resistivity from the mesh at the positions where\nthe prior data are available. To this end, we create a list of positions\n(``pg.Pos`` class) and use a forward operator that picks the values from a\nmodel vector according to the cell where the position is in. See the\nregularization tutorial for details about that.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "posVec = [pg.Pos(pos) for pos in zip(x, z)]\npara = pg.Mesh(mgr.paraDomain)  # make a copy\npara.setCellMarkers(pg.IVector(para.cellCount()))\nfopDP = PriorModelling(para, posVec)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now use it to retrieve the model values, store it and plot it along\nwith the measured values.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nax.semilogx(r, z, label=\"borehole\")\nresSmooth = fopDP(mgr.model)\nax.semilogx(resSmooth, z, label=\"ERT\")\nax.set_xlabel(r\"$\\rho$ ($\\Omega$m)\")\nax.set_ylabel(\"depth (m)\")\nax.grid(True)\nax.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As alternative to smoothness, we can use a geostatistic model. The vertical\nrange can be well estimated from the DP data using a variogram analysis, we\nguess 8m. For the horizontal one, we can only guess a ten times higher value.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mgr.inv.setRegularization(2, correlationLengths=[40, 4])\nmgr.invert()\nax, cb = mgr.showResult(**kw)\ndraw1DColumn(ax, **colkw, **kw)\nresGeo = fopDP(mgr.model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's compare the three resistivity soundings with the ground truth.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nax.semilogx(r, z, label=\"borehole\")\nax.semilogx(resSmooth, z, label=\"ERT smooth\")\n# ax.semilogx(res2, z, label=\"ERT aniso\")\nax.semilogx(resGeo, z, label=\"ERT geostat\")\nax.set_xlabel(r\"$\\rho$ ($\\Omega$m)\")\nax.set_ylabel(\"depth (m)\")\nax.grid()\nax.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The anisotropic regularization starts to see the good conductor, but only\nthe geostatistical regularization operator is able to retrieve values that\nare close to the direct push. Both show the conductor too deep.\n\nOne alternative could be to use the interfaces as structural constraints in\nthe neighborhood of the borehole. See ERT with structural constraints example\n\nAs the DP data is not only good for comparison, we want to use its values as\ndata in inversion. This is easily accomplished by taking the mapping operator\nthat we already use for interpolation as a forward operator.\n\nWe set up an inversion with this mesh, logarithmic transformations and\ninvert the model.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inv = pg.Inversion(fop=fopDP, verbose=True)\ninv.mesh = para\ntLog = pg.trans.TransLog()\ninv.modelTrans = tLog\ninv.dataTrans = tLog\ninv.setRegularization(correlationLengths=[40, 4])\nmodel = inv.run(r, relativeError=0.2)\nax, cb = pg.show(para, model, **kw)\ndraw1DColumn(ax, **colkw, **kw)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Apparently, the geostatistical operator can be used to extrapolate\nvalues with given assumptions.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Joint inversion of ERT and DP data\n\nWe now use the framework ``JointModelling`` to combine the ERT and the\nDP forward operators. So we set up a new ERT modelling operator and join\nit with ``fopDP``.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# fopERT = ert.ERTModelling()\n# fopERT.setMesh(mesh)\n# fopERT.setData(data) # not necessary as done by JointModelling\n# fopJoint = JointModelling([fopERT, fopDP])\nfopJoint = JointModelling([mgr.fop, fopDP])\n# fopJoint.setMesh(para)\nfopJoint.setData([data, pg.Vector(r)])  # needs to have .size() attribute!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first test the joint forward operator. We create a modelling vector\nof constant resistivity and distribute the model response into the two\nparts that can be looked at individually.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model = pg.Vector(para.cellCount(), 100)\nresponse = fopJoint(model)\nrespERT = response[:data.size()]\nrespDP = response[data.size():]\nprint(respDP)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The jacobian can be created and looked up by\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fopJoint.createJacobian(model)  # works\nJ = fopJoint.jacobian()\nprint(J)  # wrong type\nax, cb = pg.show(J)\nprint(J.mat(0), J.mat(1))  # ERT and DP part\nax, cb = pg.show(J.mat(1), markersize=4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the joint inversion, concatenate the data and error vectors, create a new\ninversion instance, set logarithmic transformations and run the inversion.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dataVec = np.concatenate((data[\"rhoa\"], r))\nerrorVec = np.concatenate((data[\"err\"], np.ones_like(r)*0.2))\ninv = pg.Inversion(fop=fopJoint, verbose=True)\ntransLog = pg.trans.TransLog()\ninv.modelTrans = transLog\ninv.dataTrans = transLog\ninv.run(dataVec, errorVec, startModel=model)\nax, cb = pg.show(para, inv.model, **kw)\ndraw1DColumn(ax, **colkw, **kw)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We have a local improvement of the model in the neighborhood of the\nborehole. Now we want to use geostatistics to get them further into the\nmodel.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inv.setRegularization(2, correlationLengths=[40, 4])\nmodel = inv.run(dataVec, errorVec, startModel=model)\nax, cb = pg.show(para, model, **kw)\ndraw1DColumn(ax, **colkw, **kw)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This model much better resembles the subsurface from all data and our\nexpectations to it.\n\nWe split the model response in the ERT part and the DP part. The first\nis shown as misfit.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "respERT = inv.response[:data.size()]\nmisfit = - respERT / data[\"rhoa\"] * 100 + 100\n# ax, cb = ert.show(data, misfit, cMap=\"bwr\", cMin=-5, cMax=5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The second is shown as depth profile.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "respDP = inv.response[data.size():]\nfig, ax = plt.subplots()\nax.semilogx(r, z, label=\"borehole\")\n# resMesh = pg.interpolate(srcMesh=para, inVec=inv.model, destPos=posVec)\n# ax.semilogx(resMesh, z, label=\"ERT+DP\")\nax.semilogx(respDP, z, label=\"response\")\nax.grid(True)\nax.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The model response can much better resemble the given data compared to\npure interpolation. However, at shallow depth there is some inconsistency between\nthe ERT data and the borehole data.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>Take-away messages\n\n   -  (ERT) data inversion is highly ambiguous, particularly for hidden layers\n   -  prior data can help to improve regularization\n   -  point data improve images, but only locally with smoothness constraints\n   -  geostatistical regularization can extrapolate point data\n   -  incorporation of prior data with geostatistic regularization is best</p></div>\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pg.info(\"Example finished. You can now close the window.\")"
      ]
    }
  ],
  "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.11.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}