notebooks/DivJs_inversion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file already exists, new file is called /Users/sgkang/Projects/miller2018_Tongariro_SP_inversion/notebooks/miller_et_al_2018.tar.gz\n",
      "Downloading https://storage.googleapis.com/simpeg/miller_et_al_2018.tar.gz\n",
      "   saved to: /Users/sgkang/Projects/miller2018_Tongariro_SP_inversion/notebooks/miller_et_al_2018.tar.gz\n",
      "Download completed!\n"
     ]
    }
   ],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Script to invert Tongariro SP data using divergence of streaming current density, q.\n",
    "\n",
    "from Soegi Kang, Dominique Fournier, UBC, Craig Miller, GNS\n",
    "\n",
    "https://github.com/simpeg\n",
    "\n",
    "Please use the miller_tongariro_2018 branch.  32Gb ram are recommended to run the model.\n",
    "\n",
    "\"\"\"\n",
    "from SimPEG import Utils, Maps, Mesh\n",
    "from SimPEG.EM.Static import SP\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from SimPEG.EM.Static.Utils import StaticUtils\n",
    "import os\n",
    "from sklearn.metrics import mean_squared_error as mse\n",
    "from matplotlib.ticker import AutoMinorLocator\n",
    "import colorcet as cc\n",
    "from SimPEG.Utils import download\n",
    "import tarfile\n",
    "\n",
    "tar_folder = download( \"https://storage.googleapis.com/simpeg/miller_et_al_2018.tar.gz\")  \n",
    "tar = tarfile.open(tar_folder, \"r:*\")\n",
    "tar.extractall()\n",
    "tar.close()\n",
    "\n",
    "datadir =\"./miller_et_al_2018/\"\n",
    "\n",
    "plt.close('all')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load files\n",
    "\n",
    "topo = np.loadtxt(datadir+'topoxyz.dat', skiprows=1)\n",
    "data = np.loadtxt(datadir+'Tongariro_SP_data.dat')\n",
    "\n",
    "cond_model_file = datadir+'MT_model_for_SP_inv.con'\n",
    "mesh_file = datadir+'Tongariro_SP.msh'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataInv=data\n",
    "topo = topo[(topo[:,0] < 1833000) & (topo[:,0] > 1824000)]\n",
    "topo = topo[(topo[:,1] > 5662000) & (topo[:,1] < 5670000)]\n",
    "topoInv = topo\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  ---- 3-D TensorMesh ----  \n",
      "   x0: 1824200.00\n",
      "   y0: 5662200.00\n",
      "   z0: -850.13\n",
      "  nCx: 77\n",
      "  nCy: 68\n",
      "  nCz: 86\n",
      "   hx: 207.00,  173.00,  144.00,  120.00,  70*100.00,  120.00,  144.00,  173.00,\n",
      "   hy: 207.00,  173.00,  144.00,  120.00,  60*100.00,  130.00,  144.00,  173.00,  207.00,\n",
      "   hz: 344.65,  265.11,  203.93,  156.87,  120.67,  92.82,  71.40,  54.92,  42.25,  32.50,  60*25.00,  16*20.00,\n",
      "width in X direction 8081.0 m\n",
      "width in Y direction 7298.0 m\n",
      "width in Z direction 3205.133662 m\n",
      "number of cells 450296\n"
     ]
    }
   ],
   "source": [
    "#set up the mesh\n",
    "mesh = Mesh.TensorMesh.readUBC(mesh_file)\n",
    "\n",
    "print( mesh)\n",
    "print( 'width in X direction', mesh.hx.sum(), 'm')\n",
    "print( 'width in Y direction', mesh.hy.sum(), 'm')\n",
    "print( 'width in Z direction', mesh.hz.sum(), 'm')\n",
    "print( 'number of cells', mesh.nC)\n",
    "\n",
    "xref, yref = 1830000, 5666000\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sigma is background conductivity (S/m)\n",
    "\n",
    "#read in model file\n",
    "resistivity = mesh.readModelUBC(cond_model_file) #read background conductivity model from file\n",
    "sigma = 1/resistivity #change to S/m from ohmm\n",
    "actind = sigma != 1e-8\n",
    "#sigma = np.ones(mesh.nC)*1e-3 # sets uniform background\n",
    "sigma[~actind] = 1e-8 #sigma of cells above topo = 1e-8\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set active cells\n",
    "actMap = Maps.InjectActiveCells(mesh, actind, -100.) #makes cell below topo active and sets cells above topo to 0\n",
    "mesh2D, topoCC = StaticUtils.gettopoCC(mesh, actind) #~actind = opposite of actind, ie swaps True for False\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#drape points onto topo\n",
    "xyzlocInv = StaticUtils.drapeTopotoLoc(mesh, dataInv[:,:2], actind=actind) #drape data onto topo.\n",
    "xyzlocInv = np.c_[xyzlocInv[:,:2], xyzlocInv[:,2]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAADkNJREFUeJzt3W+MXNddxvHvg9MEVCrFbtbBSmrsSFaVvGlaraJIkRAkbZImCBupQalQscDIEqKoCBC45E1ReZEgQQsSApkmwiBoEgJRrLSQGpOoQoK0Ns1f3GAnmGJsxQ5JoLwJuP3xYq7LZtn1/NmZnZ2z3480mnvP3PH8zt3xs2fPvXMnVYUkafZ917QLkCSNh4EuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJasQlq/liV1xxRW3btm01X1KSZt7Ro0dfq6q5ftutaqBv27aNI0eOrOZLStLMS/Ivg2znlIskNWKgEXqSk8A3gW8B56tqPskm4CFgG3AS+LGqemMyZUqS+hlmhP5DVXV9Vc136/uAw1W1AzjcrUuSpmQlUy47gQPd8gFg18rLkSSNatBAL+BLSY4m2du1XVlVZwC6+81LPTHJ3iRHkhw5d+7cyiuWJC1p0LNcbqqq00k2A4eSfH3QF6iq/cB+gPn5eb9NQ5ImZKARelWd7u7PAo8CNwCvJtkC0N2fnVSRkqT++gZ6kncmedeFZeBW4AXgILC722w38NikipQk9TfIlMuVwKNJLmz/p1X1V0m+CjycZA/wDeCuyZUpSeqnb6BX1SvA+5Zo/3fglkkUJa0V2/Z94TvLJ++9c4qVSP35SVFJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhqxqt9YJM2CheeeS7PEEbokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNWLgQE+yIcnXkjzerW9P8nSS40keSnLp5MqUJPUzzAj9E8CxBev3AZ+pqh3AG8CecRYmSRrOQIGe5GrgTuBz3XqAm4FHuk0OALsmUaAkaTCDjtA/C/wy8O1u/d3Am1V1vls/BVw15tokSUPoG+hJfhg4W1VHFzYvsWkt8/y9SY4kOXLu3LkRy5Qk9TPICP0m4EeSnAQepDfV8lng8iSXdNtcDZxe6slVtb+q5qtqfm5ubgwlS5KW0jfQq+qTVXV1VW0D7gb+pqp+HHgS+Ei32W7gsYlVKUnqayXnof8K8AtJTtCbU79/PCVJkkZxSf9N/k9VPQU81S2/Atww/pIkSaPwk6KS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJasRQF+eS1rNt+77wneWT9945xUqkpTlCl6RGGOiS1AgDXZIa4Ry61i3nxNUaR+iS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXC89Al3n5OujSrHKFLUiMMdElqhIEuSY0w0CWpER4UlUbghb20FjlCl6RGGOiS1AgDXZIa0TfQk3x3kq8keTbJi0l+rWvfnuTpJMeTPJTk0smXK0laziAj9LeAm6vqfcD1wO1JbgTuAz5TVTuAN4A9kytTktRP30Cvnv/qVt/R3Qq4GXikaz8A7JpIhZKkgQw0h55kQ5JngLPAIeBl4M2qOt9tcgq4ajIlSpIGMdB56FX1LeD6JJcDjwLXLrXZUs9NshfYC7B169YRy5RG5znjWi+GOsulqt4EngJuBC5PcuEXwtXA6WWes7+q5qtqfm5ubiW1SpIuYpCzXOa6kTlJvgf4IHAMeBL4SLfZbuCxSRUpSepvkCmXLcCBJBvo/QJ4uKoeT/KPwINJfh34GnD/BOuUJPXRN9Cr6jng/Uu0vwLcMImiJEnD85OiktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEb4JdHSCg168S8vEqZJc4QuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wmu5aF1ZeD0VqTWO0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmN8INF0hgt/uCSXwat1eQIXZIaYaBLUiMMdElqhHPo0gR5MTCtpr4j9CTvSfJkkmNJXkzyia59U5JDSY539xsnX64kaTmDTLmcB36xqq4FbgR+Nsl1wD7gcFXtAA5365KkKekb6FV1pqr+oVv+JnAMuArYCRzoNjsA7JpUkZKk/oY6KJpkG/B+4Gngyqo6A73QBzYv85y9SY4kOXLu3LmVVStJWtbAgZ7ke4E/B36+qv5z0OdV1f6qmq+q+bm5uVFqlCQNYKBAT/IOemH+J1X1F13zq0m2dI9vAc5OpkRJ0iAGOcslwP3Asar6rQUPHQR2d8u7gcfGX54kaVCDnId+E/Ax4Pkkz3RtvwrcCzycZA/wDeCuyZQoSRpE30Cvqr8FsszDt4y3HEnSqPzovyQ1wkCXpEYY6JLUCC/OpWYsvBCWXyyh9cgRuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjfA8dDXJL2fWeuQIXZIaYaBLUiMMdElqhHPo0pR5DRqNiyN0SWqEgS5JjTDQJakRBrokNcKDoppps/oBolmtW2ubI3RJaoSBLkmNMNAlqRHOoWvmOP8sLc0RuiQ1wkCXpEYY6JLUCOfQpRngBbw0CEfoktQIA12SGmGgS1IjDHRJakTfQE/yQJKzSV5Y0LYpyaEkx7v7jZMtU5LUzyAj9D8Ebl/Utg84XFU7gMPduiRpivoGelV9GXh9UfNO4EC3fADYNea6JElDGnUO/cqqOgPQ3W8eX0mSpFFM/INFSfYCewG2bt066ZdTo7wgl9TfqCP0V5NsAejuzy63YVXtr6r5qpqfm5sb8eUkSf2MGugHgd3d8m7gsfGUI0ka1SCnLX4e+DvgvUlOJdkD3At8KMlx4EPduiRpivrOoVfVR5d56JYx1yJJWgE/KSpJjTDQJakRBrokNcIvuNCa5bnn0nAcoUtSIwx0SWqEgS5JjXAOXZoxfmG0luMIXZIaYaBLUiMMdElqRFNz6Ksxt+j85eiW23eeb760QfaL70ct5AhdkhphoEtSIwx0SWqEgS5JjWjqoOgkrMYBu2EPbHkgrF0eINZKOEKXpEYY6JLUCANdkhrhHPoKrOW57EnUNs5/07liafwcoUtSIwx0SWqEgS5JjZiZOfRxzd+u9gW8lrMW5txXsi8W93Et9EdLm8T/nYVW+2e/lo9dLTSNOh2hS1IjDHRJaoSBLkmNmJk59ElYi+dVDztPOeyXIKxkG82O5X6eK7lu0Epedy1+4cyszMUPwxG6JDXCQJekRhjoktQIA12SGrGig6JJbgd+G9gAfK6q7h1LVUOY9MHISWntIORa+dCJBrfcQcFxHQi92DaT/iKXlfRhkH2xkpMUJmnkEXqSDcDvAh8GrgM+muS6cRUmSRrOSqZcbgBOVNUrVfXfwIPAzvGUJUka1koC/SrgXxesn+raJElTkKoa7YnJXcBtVfXT3frHgBuq6ucWbbcX2Nutvhd4afRy14wrgNemXcSU2Pf1ab32fa30+/uraq7fRis5KHoKeM+C9auB04s3qqr9wP4VvM6ak+RIVc1Pu45psO/2fT2ZtX6vZMrlq8COJNuTXArcDRwcT1mSpGGNPEKvqvNJPg48Qe+0xQeq6sWxVSZJGsqKzkOvqi8CXxxTLbOkqSmkIdn39Wm99n2m+j3yQVFJ0triR/8lqREG+hKS3JXkxSTfTjK/6LFPJjmR5KUkty1ov71rO5Fk34L27UmeTnI8yUPdAeSZkORTSf4tyTPd7Y4Fjw21H2ZZi31aLMnJJM93P+cjXdumJIe69+6hJBu79iT5nW5/PJfkA9OtfjhJHkhyNskLC9qG7muS3d32x5PsnkZf/p+q8rboBlxL75z5p4D5Be3XAc8ClwHbgZfpHRDe0C1fA1zabXNd95yHgbu75d8Hfmba/RtiP3wK+KUl2ofeD7N6a7FPy/TzJHDForbfAPZ1y/uA+7rlO4C/BALcCDw97fqH7OsPAB8AXhi1r8Am4JXufmO3vHHafXOEvoSqOlZVS30AaifwYFW9VVX/DJygdwmEJS+DkCTAzcAj3fMPALsm34OJG2o/TLHOcWixT4PaSe89C29/7+4E/qh6/h64PMmWaRQ4iqr6MvD6ouZh+3obcKiqXq+qN4BDwO2Tr/7iDPThLHe5g+Xa3w28WVXnF7XPko93f2o+cOHPUIbfD7OsxT4tpYAvJTnafbob4MqqOgPQ3W/u2lvcJ8P2dU3ug3X7naJJ/hr4viUeuqeqHlvuaUu0FUv/YqyLbL9mXGw/AL8HfJpezZ8GfhP4KYbfD7Nszf8Mx+SmqjqdZDNwKMnXL7LtetknsHxf1+Q+WLeBXlUfHOFpF7vcwVLtr9H7E+2SbpS+5OURpmnQ/ZDkD4DHu9Vh98MsG+gSF7Ouqk5392eTPEpvqunVJFuq6kw3zXC227zFfTJsX08BP7io/alVqPOinHIZzkHg7iSXJdkO7AC+wjKXQaje0ZMngY90z98NLDf6X3MWzYv+KHDhrICh9sNq1jwBLfbpbZK8M8m7LiwDt9L7WR+k956Ft793DwI/0Z0BciPwHxemK2bYsH19Arg1ycZuKvLWrm26pn1Udi3e6IXXKeAt4FXgiQWP3UPvrIeXgA8vaL8D+KfusXsWtF9DL+xOAH8GXDbt/g2xH/4YeB54jt4be8uo+2GWby32aVH/rqF39s6zwIsX+kjvGNBh4Hh3v6lrD70vt3m5e3/MT6v2Efv7eeAM8D/d//M9o/SV3vTjie72k9PuV1X5SVFJaoVTLpLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RG/C8GsPpPX33tNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113338b70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#  Inversion runs in the next sections.\n",
    "# Setup the problem\n",
    "\n",
    "from pymatsolver import PardisoSolver\n",
    "\n",
    "ind = np.argmin(abs(xyzlocInv[:, 0]-xref)+abs(xyzlocInv[:, 1]-yref))\n",
    "inds = np.ones(xyzlocInv.shape[0], dtype='bool')\n",
    "inds[ind] = False\n",
    "xyzM = xyzlocInv[inds, :]\n",
    "xyzN = np.atleast_2d(xyzlocInv[ind]).repeat(xyzM.shape[0], axis=0)\n",
    "#wires = Maps.Wires(('jsx', actMap.nP), ('jsy', actMap.nP), ('jsz', actMap.nP))\n",
    "prob = SP.Problem_CC_Jstore(mesh, sigma=sigma, qMap=actMap, Solver=PardisoSolver)\n",
    "rx = SP.Rx.Dipole(xyzM, xyzN)\n",
    "src = SP.Src.StreamingCurrents([rx], L=np.ones(mesh.nC), mesh=mesh, modelType=\"CurrentSource\")\n",
    "survey = SP.SurveySP_store([src])\n",
    "survey.pair(prob)\n",
    "\n",
    "survey.dobs = (dataInv[1:,2]-dataInv[0,2])  # mV\n",
    "\n",
    "\n",
    "plt.hist(survey.dobs, bins=100)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from SimPEG import (Mesh, Maps, DataMisfit, Regularization,\n",
    "                    Optimization, Inversion, InvProblem, Directives)\n",
    "\n",
    "survey.std = 30 #survey std dev in mV\n",
    "survey.eps = 60 # survey noise floor in mV\n",
    "\n",
    "\n",
    "dmisfit = DataMisfit.l2_DataMisfit(survey)\n",
    "regmap = Maps.IdentityMap(nP = actMap.nP)\n",
    "reg = Regularization.Sparse(mesh, mapping=regmap, indActive=actind)\n",
    "\n",
    "reg.alpha_s = 1\n",
    "\n",
    "reg.norms = [0, 1, 1, 1]\n",
    "\n",
    "opt = Optimization.ProjectedGNCG(maxIter=100, maxIterCG=10, tolX=1e-20, tolF=1e-20)\n",
    "\n",
    "opt.maxIterLS = 20\n",
    "\n",
    "invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)\n",
    "target = Directives.TargetMisfit()\n",
    "\n",
    "# Create an inversion object\n",
    "update_Jacobi = Directives.UpdatePreconditioner()\n",
    "update_SensWeight = Directives.UpdateSensitivityWeights()\n",
    "\n",
    "IRLS = Directives.Update_IRLS(\n",
    "        f_min_change=1e-3, minGNiter=2, maxIRLSiter = 10,\n",
    "        coolingFactor=8, coolingRate=3\n",
    "        )\n",
    "\n",
    "\n",
    "betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1.)\n",
    "updateprecond = Directives.UpdatePreconditioner()\n",
    "\n",
    "#LP only\n",
    "inv = Inversion.BaseInversion(invProb, directiveList=[betaest, IRLS, updateprecond, update_SensWeight])\n",
    "\n",
    "#save iterations\n",
    "prob.counter = opt.counter = Utils.Counter()\n",
    "opt.LSshorten = 0.5\n",
    "opt.remember('xc')\n",
    "\n",
    "#creates an initial starting model m0\n",
    "m0 = np.ones(actMap.nP)*0.\n",
    "reg.mref = m0\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
      "    ***Done using same Solver and solverOpts as the problem***\n"
     ]
    }
   ],
   "source": [
    "#run the inversion\n",
    "mopt_Lp = inv.run(m0)\n",
    "Lp_amp = actMap*mopt_Lp\n",
    "\n",
    "#recall a particular iteration, comment out if not needed.\n",
    "xc = opt.recall('xc')\n",
    "mopt_Lp = xc[8]\n",
    "Lp_amp = actMap*mopt_Lp\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save models and calc data\n",
    "\n",
    "mesh.writeModelUBC('Tongairo_C_Lp_amp.mod', Lp_amp)\n",
    "\n",
    "#Write calculated data (Lp model)\n",
    "Lp_calculated = survey.dpred(mopt_Lp)\n",
    "np.savetxt('Lp_calculated_data.dat', Lp_calculated, fmt='%.4f')\n",
    "\n",
    "\n",
    "RMS = np.sqrt(mse(Lp_calculated, survey.dobs))\n",
    "print('RMS = %.3f mV' %RMS)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot some L2 model obs and calc\n",
    "plt.figure()\n",
    "plt.plot(survey.dobs, '.', label = 'obs')\n",
    "plt.plot(invProb.dpred, label = 'calc')\n",
    "plt.legend()\n",
    "plt.ylabel('mV')\n",
    "plt.xlabel('data point index')\n",
    "plt.annotate('B', xy=(0, 1000), fontsize=14)\n",
    "plt.title('$q$ observed vs calculated.  RMS = %.0f mV' %RMS)\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(mopt_Lp, bins=100, label='Lp model')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract core mesh\n",
    "temp = Lp_amp.copy()\n",
    "temp[~actind] = np.nan\n",
    "vmin, vmax = temp[actind].min(), temp[actind].max()\n",
    "\n",
    "xmin, xmax = dataInv[:,0].min(), dataInv[:,0].max()\n",
    "ymin, ymax = dataInv[:,1].min(), dataInv[:,1].max()\n",
    "zmin, zmax = 700, topo[:,2].max()\n",
    "xyzlim = np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])\n",
    "ind_core, mesh_core = Utils.ExtractCoreMesh(xyzlim, mesh)\n",
    "\n",
    "ja_est_core = Utils.mkvc(temp[ind_core])\n",
    "actind_core = actind[ind_core]\n",
    "jamp_min, jamp_max = ja_est_core[:mesh_core.nC][actind_core].min(), ja_est_core[:mesh_core.nC][actind_core].max()\n",
    "jamp_lp_core = temp[ind_core]\n",
    "jamp_lp_core[~actind_core] = np.nan\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot map views\n",
    "fig, ((ax0,ax1),(ax2,ax3)) = plt.subplots(nrows=2, ncols=2, figsize = (10, 10))\n",
    "\n",
    "fig.subplots_adjust(wspace=0, hspace=0)\n",
    "axs = [ax0, ax1, ax2, ax3]\n",
    "layers = [1500, 1200, 1000, 800]\n",
    "\n",
    "\n",
    "clim=(-0.00005, 0.00005)\n",
    "for i_layer, ax in enumerate(axs):   \n",
    "    jamp_lp_core[~actind_core] = np.nan\n",
    "    indz = int(np.argmin(abs(mesh_core.vectorCCz-layers[i_layer])))\n",
    "    cbar = mesh_core.plotSlice(jamp_lp_core, ind=indz, normal='Z', clim=clim, ax=ax)\n",
    "    \n",
    "    minorLocator = AutoMinorLocator(4)\n",
    "    ax.xaxis.set_minor_locator(minorLocator)\n",
    "    ax.yaxis.set_minor_locator(minorLocator)\n",
    "\n",
    "    ax.tick_params(which='both', width=2, direction='in')\n",
    "    ax.tick_params(which='major', length=10, direction='in', top= True, right=True)\n",
    "    ax.tick_params(which='minor', length=7, direction='in', top = True)\n",
    "    \n",
    "    ax.set_title('')\n",
    "    ax.set_xlabel('')\n",
    "    ax.set_ylabel('')\n",
    "    ax.set_aspect('equal')\n",
    "    ax.set_ylim(5663000, 5669000)\n",
    "    ax.set_xlim(1825094, 1831794)    \n",
    "    ax.set_title((\"%.0f\")%(layers[i_layer]))\n",
    "    ax.set_xticks(np.arange(1826000, 1832000, 2000))\n",
    "    ax.set_yticks(np.arange(5663000, 5670000, 2000))\n",
    "    \n",
    "    if ~np.logical_or(i_layer==2,i_layer==3):\n",
    "        ax.set_xticklabels([])\n",
    "    if np.logical_or(i_layer==1,i_layer==3):\n",
    "        ax.set_yticklabels({})\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot section views\n",
    "fig, ((ax0,ax1),(ax2,ax3),(ax4,ax5)) = plt.subplots(nrows=3, ncols=2, figsize = (10, 10))\n",
    "\n",
    "fig.delaxes(ax5)\n",
    "\n",
    "fig.subplots_adjust(wspace=0, hspace=0)\n",
    "axs = [ax0, ax1, ax2, ax3, ax4, ax5]\n",
    "layers = [5667994, 5666994, 5665994, 5664994, 5663994]\n",
    "\n",
    "stream_color = 'k'\n",
    "clim = (-0.00005, 0.00005)\n",
    "for i_layer, ax in enumerate(axs[:-1]):\n",
    "    jamp_lp_core[~actind_core] = np.nan\n",
    "    indy = int(np.argmin(abs(mesh_core.vectorCCy-layers[i_layer])))\n",
    "    cbar = mesh_core.plotSlice(jamp_lp_core, ind=indy, normal='Y', clim=clim, ax=ax)\n",
    "\n",
    "    minorLocator = AutoMinorLocator(4)\n",
    "    ax.xaxis.set_minor_locator(minorLocator)\n",
    "    ax.yaxis.set_minor_locator(minorLocator)\n",
    "\n",
    "    ax.tick_params(which='both', width=2, direction='in')\n",
    "    ax.tick_params(which='major', length=10, direction='in', top= True, right=True)\n",
    "    ax.tick_params(which='minor', length=7, direction='in', top = True)\n",
    "    \n",
    "    ax.set_title('')\n",
    "    ax.set_xlabel('')\n",
    "    ax.set_ylabel('')\n",
    "    ax.set_aspect('equal')\n",
    "    ax.set_ylim(700, 2000)\n",
    "    ax.set_xlim(1825094, 1831794)    \n",
    "    ax.set_title((\"%.0f\")%(layers[i_layer]))\n",
    "    ax.set_xticks([1826000, 1828000, 1830000])\n",
    "    ax.set_yticks([1000, 1500, 2000])\n",
    "      \n",
    "    if ~np.logical_or(i_layer==4,i_layer==3) :\n",
    "        ax.set_xticklabels([])\n",
    "    if np.logical_or(i_layer==1,i_layer==3):\n",
    "        ax.set_yticklabels({})\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write core of the L2 model to file in spherical formulation\n",
    "mesh_core.writeUBC('mesh_core.msh')\n",
    "\n",
    "jamp_lp_core[~actind_core] = -100\n",
    "mesh_core.writeModelUBC('Tongairo_C_Lp_amp_core.mod', jamp_lp_core)"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}