ERT/3dtopo/slagdump3d.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3D ERT inversion with topography\n",
    "\n",
    "We demonstrate 3D inversion with topography by using an example data set that is described in the original BERT paper (Günther et al., 2006). The data set images an industrial waste dump with 13 cross-profiles of 2m electrode distance and a long profile along the ridge with 1m electrode distance, a total of more than 5000 data with almost 600 electrodes.\n",
    "\n",
    "Günther, T., Rücker, C. & Spitzer, K. (2006): Three-dimensional modeling and inversion of dc resistivity data incorporating topography – II: Inversion. Geophys. J. Int. 166, 506-517,doi:10.1111/j.1365-246X.2006.03011.x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pygimli as pg\n",
    "import pygimli.meshtools as mt\n",
    "from pygimli.physics import ert\n",
    "from pygimli.viewer import pv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data: Sensors: 577 data: 4245, nonzero entries: ['a', 'b', 'm', 'n', 'r', 'valid']\n"
     ]
    }
   ],
   "source": [
    "data = ert.load(\"slagdump3d.ohm\")\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mesh: Nodes: 2450 Cells: 0 Boundaries: 4875\n"
     ]
    }
   ],
   "source": [
    "plc = mt.createParaMeshPLC3D(data, paraDepth=30)\n",
    "print(plc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mesh: Nodes: 5466 Cells: 24598 Boundaries: 52003\n"
     ]
    },
    {
     "data": {
      "application/javascript": "(function(root) {\n  function now() {\n    return new Date();\n  }\n\n  var force = true;\n\n  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n    root._bokeh_onload_callbacks = [];\n    root._bokeh_is_loading = undefined;\n  }\n\n  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n    root._bokeh_timeout = Date.now() + 5000;\n    root._bokeh_failed_load = false;\n  }\n\n  function run_callbacks() {\n    try {\n      root._bokeh_onload_callbacks.forEach(function(callback) {\n        if (callback != null)\n          callback();\n      });\n    } finally {\n      delete root._bokeh_onload_callbacks\n    }\n    console.debug(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(css_urls, js_urls, js_modules, callback) {\n    if (css_urls == null) css_urls = [];\n    if (js_urls == null) js_urls = [];\n    if (js_modules == null) js_modules = [];\n\n    root._bokeh_onload_callbacks.push(callback);\n    if (root._bokeh_is_loading > 0) {\n      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls.length === 0 && js_modules.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n\n    function on_load() {\n      root._bokeh_is_loading--;\n      if (root._bokeh_is_loading === 0) {\n        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n        run_callbacks()\n      }\n    }\n\n    function on_error() {\n      console.error(\"failed to load \" + url);\n    }\n\n    for (var i = 0; i < css_urls.length; i++) {\n      var url = css_urls[i];\n      const element = document.createElement(\"link\");\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.rel = \"stylesheet\";\n      element.type = \"text/css\";\n      element.href = url;\n      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n      document.body.appendChild(element);\n    }\n\n    var skip = [];\n    if (window.requirejs) {\n      window.requirejs.config({'packages': {}, 'paths': {'vtk': 'https://cdn.jsdelivr.net/npm/vtk.js@20.0.1/vtk', 'gridstack': 'https://cdn.jsdelivr.net/npm/gridstack@4.2.5/dist/gridstack-h5', 'notyf': 'https://cdn.jsdelivr.net/npm/notyf@3/notyf.min'}, 'shim': {'vtk': {'exports': 'vtk'}, 'gridstack': {'exports': 'GridStack'}}});\n      require([\"vtk\"], function() {\n\ton_load()\n      })\n      require([\"gridstack\"], function(GridStack) {\n\twindow.GridStack = GridStack\n\ton_load()\n      })\n      require([\"notyf\"], function() {\n\ton_load()\n      })\n      root._bokeh_is_loading = css_urls.length + 3;\n    } else {\n      root._bokeh_is_loading = css_urls.length + js_urls.length + js_modules.length;\n    }    if (((window['vtk'] !== undefined) && (!(window['vtk'] instanceof HTMLElement))) || window.requirejs) {\n      var urls = ['https://cdn.holoviz.org/panel/0.14.4/dist/bundled/abstractvtkplot/vtk.js@20.0.1/vtk.js'];\n      for (var i = 0; i < urls.length; i++) {\n        skip.push(urls[i])\n      }\n    }    if (((window['GridStack'] !== undefined) && (!(window['GridStack'] instanceof HTMLElement))) || window.requirejs) {\n      var urls = ['https://cdn.holoviz.org/panel/0.14.4/dist/bundled/gridstack/gridstack@4.2.5/dist/gridstack-h5.js'];\n      for (var i = 0; i < urls.length; i++) {\n        skip.push(urls[i])\n      }\n    }    if (((window['Notyf'] !== undefined) && (!(window['Notyf'] instanceof HTMLElement))) || window.requirejs) {\n      var urls = ['https://cdn.holoviz.org/panel/0.14.4/dist/bundled/notificationarea/notyf@3/notyf.min.js'];\n      for (var i = 0; i < urls.length; i++) {\n        skip.push(urls[i])\n      }\n    }    for (var i = 0; i < js_urls.length; i++) {\n      var url = js_urls[i];\n      if (skip.indexOf(url) >= 0) {\n\tif (!window.requirejs) {\n\t  on_load();\n\t}\n\tcontinue;\n      }\n      var element = document.createElement('script');\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.async = false;\n      element.src = url;\n      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.head.appendChild(element);\n    }\n    for (var i = 0; i < js_modules.length; i++) {\n      var url = js_modules[i];\n      if (skip.indexOf(url) >= 0) {\n\tif (!window.requirejs) {\n\t  on_load();\n\t}\n\tcontinue;\n      }\n      var element = document.createElement('script');\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.async = false;\n      element.src = url;\n      element.type = \"module\";\n      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.head.appendChild(element);\n    }\n    if (!js_urls.length && !js_modules.length) {\n      on_load()\n    }\n  };\n\n  function inject_raw_css(css) {\n    const element = document.createElement(\"style\");\n    element.appendChild(document.createTextNode(css));\n    document.body.appendChild(element);\n  }\n\n  var js_urls = [\"https://cdn.holoviz.org/panel/0.14.4/dist/bundled/abstractvtkplot/vtk.js@20.0.1/vtk.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-2.4.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-2.4.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-2.4.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-2.4.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-2.4.3.min.js\", \"https://unpkg.com/@holoviz/panel@0.14.4/dist/panel.min.js\"];\n  var js_modules = [];\n  var css_urls = [\"https://cdn.holoviz.org/panel/0.14.4/dist/css/alerts.css\", \"https://cdn.holoviz.org/panel/0.14.4/dist/css/card.css\", \"https://cdn.holoviz.org/panel/0.14.4/dist/css/dataframe.css\", \"https://cdn.holoviz.org/panel/0.14.4/dist/css/debugger.css\", \"https://cdn.holoviz.org/panel/0.14.4/dist/css/json.css\", \"https://cdn.holoviz.org/panel/0.14.4/dist/css/loading.css\", \"https://cdn.holoviz.org/panel/0.14.4/dist/css/markdown.css\", \"https://cdn.holoviz.org/panel/0.14.4/dist/css/widgets.css\"];\n  var inline_js = [    function(Bokeh) {\n      inject_raw_css(\"\\n    .bk.pn-loading.arc:before {\\n      background-image: url(\\\"data:image/svg+xml;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHN0eWxlPSJtYXJnaW46IGF1dG87IGJhY2tncm91bmQ6IG5vbmU7IGRpc3BsYXk6IGJsb2NrOyBzaGFwZS1yZW5kZXJpbmc6IGF1dG87IiB2aWV3Qm94PSIwIDAgMTAwIDEwMCIgcHJlc2VydmVBc3BlY3RSYXRpbz0ieE1pZFlNaWQiPiAgPGNpcmNsZSBjeD0iNTAiIGN5PSI1MCIgZmlsbD0ibm9uZSIgc3Ryb2tlPSIjYzNjM2MzIiBzdHJva2Utd2lkdGg9IjEwIiByPSIzNSIgc3Ryb2tlLWRhc2hhcnJheT0iMTY0LjkzMzYxNDMxMzQ2NDE1IDU2Ljk3Nzg3MTQzNzgyMTM4Ij4gICAgPGFuaW1hdGVUcmFuc2Zvcm0gYXR0cmlidXRlTmFtZT0idHJhbnNmb3JtIiB0eXBlPSJyb3RhdGUiIHJlcGVhdENvdW50PSJpbmRlZmluaXRlIiBkdXI9IjFzIiB2YWx1ZXM9IjAgNTAgNTA7MzYwIDUwIDUwIiBrZXlUaW1lcz0iMDsxIj48L2FuaW1hdGVUcmFuc2Zvcm0+ICA8L2NpcmNsZT48L3N2Zz4=\\\");\\n      background-size: auto calc(min(50%, 400px));\\n    }\\n    \");\n    },    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\nfunction(Bokeh) {} // ensure no trailing comma for IE\n  ];\n\n  function run_inline_js() {\n    if ((root.Bokeh !== undefined) || (force === true)) {\n      for (var i = 0; i < inline_js.length; i++) {\n        inline_js[i].call(root, root.Bokeh);\n      }} else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(run_inline_js, 100);\n    } else if (!root._bokeh_failed_load) {\n      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n      root._bokeh_failed_load = true;\n    }\n  }\n\n  if (root._bokeh_is_loading === 0) {\n    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(css_urls, js_urls, js_modules, function() {\n      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(window));",
      "application/vnd.holoviews_load.v0+json": ""
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": "\nif ((window.PyViz === undefined) || (window.PyViz instanceof HTMLElement)) {\n  window.PyViz = {comms: {}, comm_status:{}, kernels:{}, receivers: {}, plot_index: []}\n}\n\n\n    function JupyterCommManager() {\n    }\n\n    JupyterCommManager.prototype.register_target = function(plot_id, comm_id, msg_handler) {\n      if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n        var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n        comm_manager.register_target(comm_id, function(comm) {\n          comm.on_msg(msg_handler);\n        });\n      } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n        window.PyViz.kernels[plot_id].registerCommTarget(comm_id, function(comm) {\n          comm.onMsg = msg_handler;\n        });\n      } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n        google.colab.kernel.comms.registerTarget(comm_id, (comm) => {\n          var messages = comm.messages[Symbol.asyncIterator]();\n          function processIteratorResult(result) {\n            var message = result.value;\n            console.log(message)\n            var content = {data: message.data, comm_id};\n            var buffers = []\n            for (var buffer of message.buffers || []) {\n              buffers.push(new DataView(buffer))\n            }\n            var metadata = message.metadata || {};\n            var msg = {content, buffers, metadata}\n            msg_handler(msg);\n            return messages.next().then(processIteratorResult);\n          }\n          return messages.next().then(processIteratorResult);\n        })\n      }\n    }\n\n    JupyterCommManager.prototype.get_client_comm = function(plot_id, comm_id, msg_handler) {\n      if (comm_id in window.PyViz.comms) {\n        return window.PyViz.comms[comm_id];\n      } else if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n        var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n        var comm = comm_manager.new_comm(comm_id, {}, {}, {}, comm_id);\n        if (msg_handler) {\n          comm.on_msg(msg_handler);\n        }\n      } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n        var comm = window.PyViz.kernels[plot_id].connectToComm(comm_id);\n        comm.open();\n        if (msg_handler) {\n          comm.onMsg = msg_handler;\n        }\n      } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n        var comm_promise = google.colab.kernel.comms.open(comm_id)\n        comm_promise.then((comm) => {\n          window.PyViz.comms[comm_id] = comm;\n          if (msg_handler) {\n            var messages = comm.messages[Symbol.asyncIterator]();\n            function processIteratorResult(result) {\n              var message = result.value;\n              var content = {data: message.data};\n              var metadata = message.metadata || {comm_id};\n              var msg = {content, metadata}\n              msg_handler(msg);\n              return messages.next().then(processIteratorResult);\n            }\n            return messages.next().then(processIteratorResult);\n          }\n        }) \n        var sendClosure = (data, metadata, buffers, disposeOnDone) => {\n          return comm_promise.then((comm) => {\n            comm.send(data, metadata, buffers, disposeOnDone);\n          });\n        };\n        var comm = {\n          send: sendClosure\n        };\n      }\n      window.PyViz.comms[comm_id] = comm;\n      return comm;\n    }\n    window.PyViz.comm_manager = new JupyterCommManager();\n    \n\n\nvar JS_MIME_TYPE = 'application/javascript';\nvar HTML_MIME_TYPE = 'text/html';\nvar EXEC_MIME_TYPE = 'application/vnd.holoviews_exec.v0+json';\nvar CLASS_NAME = 'output';\n\n/**\n * Render data to the DOM node\n */\nfunction render(props, node) {\n  var div = document.createElement(\"div\");\n  var script = document.createElement(\"script\");\n  node.appendChild(div);\n  node.appendChild(script);\n}\n\n/**\n * Handle when a new output is added\n */\nfunction handle_add_output(event, handle) {\n  var output_area = handle.output_area;\n  var output = handle.output;\n  if ((output.data == undefined) || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n    return\n  }\n  var id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n  var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n  if (id !== undefined) {\n    var nchildren = toinsert.length;\n    var html_node = toinsert[nchildren-1].children[0];\n    html_node.innerHTML = output.data[HTML_MIME_TYPE];\n    var scripts = [];\n    var nodelist = html_node.querySelectorAll(\"script\");\n    for (var i in nodelist) {\n      if (nodelist.hasOwnProperty(i)) {\n        scripts.push(nodelist[i])\n      }\n    }\n\n    scripts.forEach( function (oldScript) {\n      var newScript = document.createElement(\"script\");\n      var attrs = [];\n      var nodemap = oldScript.attributes;\n      for (var j in nodemap) {\n        if (nodemap.hasOwnProperty(j)) {\n          attrs.push(nodemap[j])\n        }\n      }\n      attrs.forEach(function(attr) { newScript.setAttribute(attr.name, attr.value) });\n      newScript.appendChild(document.createTextNode(oldScript.innerHTML));\n      oldScript.parentNode.replaceChild(newScript, oldScript);\n    });\n    if (JS_MIME_TYPE in output.data) {\n      toinsert[nchildren-1].children[1].textContent = output.data[JS_MIME_TYPE];\n    }\n    output_area._hv_plot_id = id;\n    if ((window.Bokeh !== undefined) && (id in Bokeh.index)) {\n      window.PyViz.plot_index[id] = Bokeh.index[id];\n    } else {\n      window.PyViz.plot_index[id] = null;\n    }\n  } else if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n    var bk_div = document.createElement(\"div\");\n    bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n    var script_attrs = bk_div.children[0].attributes;\n    for (var i = 0; i < script_attrs.length; i++) {\n      toinsert[toinsert.length - 1].childNodes[1].setAttribute(script_attrs[i].name, script_attrs[i].value);\n    }\n    // store reference to server id on output_area\n    output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n  }\n}\n\n/**\n * Handle when an output is cleared or removed\n */\nfunction handle_clear_output(event, handle) {\n  var id = handle.cell.output_area._hv_plot_id;\n  var server_id = handle.cell.output_area._bokeh_server_id;\n  if (((id === undefined) || !(id in PyViz.plot_index)) && (server_id !== undefined)) { return; }\n  var comm = window.PyViz.comm_manager.get_client_comm(\"hv-extension-comm\", \"hv-extension-comm\", function () {});\n  if (server_id !== null) {\n    comm.send({event_type: 'server_delete', 'id': server_id});\n    return;\n  } else if (comm !== null) {\n    comm.send({event_type: 'delete', 'id': id});\n  }\n  delete PyViz.plot_index[id];\n  if ((window.Bokeh !== undefined) & (id in window.Bokeh.index)) {\n    var doc = window.Bokeh.index[id].model.document\n    doc.clear();\n    const i = window.Bokeh.documents.indexOf(doc);\n    if (i > -1) {\n      window.Bokeh.documents.splice(i, 1);\n    }\n  }\n}\n\n/**\n * Handle kernel restart event\n */\nfunction handle_kernel_cleanup(event, handle) {\n  delete PyViz.comms[\"hv-extension-comm\"];\n  window.PyViz.plot_index = {}\n}\n\n/**\n * Handle update_display_data messages\n */\nfunction handle_update_output(event, handle) {\n  handle_clear_output(event, {cell: {output_area: handle.output_area}})\n  handle_add_output(event, handle)\n}\n\nfunction register_renderer(events, OutputArea) {\n  function append_mime(data, metadata, element) {\n    // create a DOM node to render to\n    var toinsert = this.create_output_subarea(\n    metadata,\n    CLASS_NAME,\n    EXEC_MIME_TYPE\n    );\n    this.keyboard_manager.register_events(toinsert);\n    // Render to node\n    var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n    render(props, toinsert[0]);\n    element.append(toinsert);\n    return toinsert\n  }\n\n  events.on('output_added.OutputArea', handle_add_output);\n  events.on('output_updated.OutputArea', handle_update_output);\n  events.on('clear_output.CodeCell', handle_clear_output);\n  events.on('delete.Cell', handle_clear_output);\n  events.on('kernel_ready.Kernel', handle_kernel_cleanup);\n\n  OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n    safe: true,\n    index: 0\n  });\n}\n\nif (window.Jupyter !== undefined) {\n  try {\n    var events = require('base/js/events');\n    var OutputArea = require('notebook/js/outputarea').OutputArea;\n    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n      register_renderer(events, OutputArea);\n    }\n  } catch(err) {\n  }\n}\n",
      "application/vnd.holoviews_load.v0+json": ""
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<style>.bk-root, .bk-root .bk:before, .bk-root .bk:after {\n",
       "  font-family: var(--jp-ui-font-size1);\n",
       "  font-size: var(--jp-ui-font-size1);\n",
       "  color: var(--jp-ui-font-color1);\n",
       "}\n",
       "</style>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bc41d2fe59684501939675b0b4e98a39",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "BokehModel(combine_events=True, render_bundle={'docs_json': {'d420efc6-4c9f-4e61-a0fc-1b093d0a3e4b': {'defs': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(<pyvista.plotting.plotting.Plotter at 0x2013376ae80>, None)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mesh = mt.createMesh(plc, quality=1.3)\n",
    "print(mesh)\n",
    "pg.show(mesh, notebook=True)  # backend=\"jupyter\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "28/08/23 - 11:14:08 - pyGIMLi - \u001b[0;33;49mWARNING\u001b[0m - Sensor z-coordinates not equal. Is there topography?\n",
      "28/08/23 - 11:14:08 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - mesh Mesh: Nodes: 5466 Cells: 24598 Boundaries: 52003\n",
      "28/08/23 - 11:14:09 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - p2 refine Mesh: Nodes: 38336 Cells: 24598 Boundaries: 52003\n",
      "28/08/23 - 11:14:09 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Calculate numerical geometric factors.\n",
      "28/08/23 - 11:14:22 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Save RVector binary\n",
      "28/08/23 - 11:14:22 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Cache stored: C:\\Guenther.T\\Anwendungsdaten\\pygimli\\Cache\\13469656690945789773\n"
     ]
    }
   ],
   "source": [
    "data[\"k\"] = ert.createGeometricFactors(data, mesh=mesh, h2=False, p2=True)\n",
    "# k = ert.createGeometricFactors(data, numerical=True)  # fails\n",
    "data[\"rhoa\"] = data[\"r\"] * data[\"k\"]\n",
    "# pg.plt.hist(pg.log10(data[\"rhoa\"]), bins=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Program Files\\software\\Anaconda\\envs\\pg14\\lib\\site-packages\\numpy\\core\\fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "c:\\Program Files\\software\\Anaconda\\envs\\pg14\\lib\\site-packages\\numpy\\core\\_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "28/08/23 - 11:14:22 - pyGIMLi - \u001b[0;33;49mWARNING\u001b[0m - Something gone wrong while drawing data. Try fallback with equidistant electrodes.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data: Sensors: 38 data: 222, nonzero entries: ['a', 'b', 'k', 'm', 'n', 'r', 'rhoa', 'valid']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<Axes: >, <matplotlib.colorbar.Colorbar at 0x2014e972c40>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAG0CAYAAAAmZLNuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABF6ElEQVR4nO3de3hU1b3/8c8kmSQQkhGEEEJCIIkJKooXeioqTIKIoV6woAUVC5Xi+VFAjq1tRT1HesVesLa1lWoTlJajYsUqiEcsZjAQIkEI4I2ggqCgWAi5kQuTrN8fkC2T7MllJyRg3q/nmSfstfd3rbUvyXzYmZm4jDFGAAAACBDS1RMAAAA4HRGSAAAAbBCSAAAAbBCSAAAAbBCSAAAAbBCSAAAAbBCSAAAAbIR19QQ6U319vfbv36/o6Gi5XK6ung4AAGgFY4zKy8sVHx+vkJDOu7/TrULS/v37lZiY2NXTAAAADuzbt08JCQmdNl63CknR0dGSpM3XZ6iXu2273n9SnKMxj7yy31Fd7x9nOqqTJPX0OCrbE+7sw9cHV9Y4qvv9v89yVHdXxTZHdZIUct7FjuoO3/e/juoKX3N2x/LQF35HdXEJbkd1g5Kd/yhIWTzGWWGPsxyV7QmpclRX9O8ejupuDD/kqO7Hn/Z1VCdJvx7q7Lop8lc7qjvnV885qnOn9HZUF3bhQEd1xzZ/7KhOkjbffK2justdzq43uXs6Kqv19HFUZ1TvqK5yzuOO6iSpttrZHZ21K2rbXFMlv+5RvvU83lm6VUhq+BVbL3eYot1tezKJ6RnuaMz6cGeHOCY60lGdJKmnsyeD6AhnISnG4a3PyGpnP0RiTISjOkkKiXE2pt/heewZ4uzJ7qjD3wZHhTgLSb3CnP8oiIl2dr05vk4d3mnvWePweguvdFQXEeVsPEmKiXF2AfQ65qwuxuH17e7h7HoLi3L2PXws0tl4khTl8Dp1eCqkcGfj1Tr8GeU0JIU4PPeSVFsX6qiuh8vBXE88PXX2S2V44TYAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAICNdoWkxYsXKzo6Wn7/l5/pUlFRIbfbrVGjRgVsm5eXJ5fLpeLiYqstPz9foaGhysrKatV4K1as0DXXXKO+ffvK5XKpqKioPdMHAAAIql0hKTMzUxUVFdq8ebPVlpeXp7i4OBUWFuro0aNWu8/nU3x8vNLS0qy2nJwczZ07V+vXr9fevXtbHK+yslJXXHGFHnroofZMGwAAoEXtCknp6emKj4+Xz+ez2nw+nyZMmKCUlBTl5+cHtGdmfvkp0pWVlVq+fLlmzZql6667Tk8++WSL491+++36n//5H40dO7Y90wYAAGhRu1+TlJGRodzcXGs5NzdXGRkZ8nq9Vnttba02btwYEJKeffZZpaenKz09XVOnTtWSJUtkjLNPfA6mpqZGZWVlAQ8AAIDW6JCQtGHDBvn9fpWXl2vr1q0aPXq0vF6vdYepoKBAVVVVASEpOztbU6dOlSRlZWWpoqJCa9eube90AixcuFAej8d68MdtAQBAa7U7JGVmZqqyslKFhYXKy8tTWlqaYmNj5fV6VVhYqMrKSvl8Pg0aNEjJycmSpJ07d2rTpk2aMmWKJCksLEyTJ09WTk5Oe6cTYP78+SotLbUe+/bt69D+AQDAV1e7/8BtamqqEhISlJubq5KSEnm9XklSXFychgwZog0bNig3N1djxnz5l8Kzs7Pl9/s1cOCXfwnaGCO3262SkhL17u3sL0s3FhERoYgI538MFQAAdF8d8jlJmZmZ8vl88vl8ysjIsNq9Xq9effVVFRQUWL9q8/v9Wrp0qRYtWqSioiLrsW3bNiUlJWnZsmUdMSUAAIB2afedJOl4SJo9e7aOHTtm3UmSjoekWbNmqbq62gpJq1atUklJiWbMmCGPxxPQz0033aTs7GzNmTPHdpzDhw9r79692r9/v6Tjv7aTjt+1iouL64hdAQAAkNSBd5KqqqqUmpqq/v37W+1er1fl5eVKSUmxXjSdnZ2tsWPHNglIkjRp0iQVFRVpy5YttuO89NJLuvjii3XttddKkqZMmaKLL75Yixcv7ojdAAAAsHTInaTBgwfbvn0/ISGhSfvKlSuD9nPJJZc0+zEA06dP1/Tp0x3PEwAAoLX4220AAAA2CEkAAAA2CEkAAAA2CEkAAAA2CEkAAAA2XKaj/6rsaaysrEwej0elpaWKiYnpsH6/qF7iqC4ytJfjMR9/r9ZR3aFql6O6b6VUO6oLczm7vEId1oWFOL+cz4lOdlRnyg44G/Dzj53V+euc1Z11lrO6dnDFneuw0Nn/394+ssdR3R/f6eGobvvrfR3VjRz3b0d1knRFf2fX+ITBzuZaZ/yO6rpC+JEvHNV93uOYozpPeKyjOr9x9vN7zT5ndd8sXu+oTpJcFzj7HnYNvMvxmHZO1fN3S7iTBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYIOQBAAAYMNljDFdPYnOUlZWJo/Ho9LSUsXExHT1dJr15M5ljmvfKXE5qsuMr3M8phPD+lQ5qvvbrl6Ox7ygd72jun49/I7qzu/d01Fd+bESR3VH/c7+3/PGAWfzlKRwh//V6h3h7FycFeHsOv38qNtR3eo9oY7q9q6KclQnSTGXO7ve/t9l5Y7qUmJqHNX1cjs7NtsOhTmq+/26sxzVSZLb7ex6mzC80lHd+EHOfr7V1Dmb58CoIY7qQl1XO6rrbF31/M2dJAAAABuEJAAAABuEJAAAABuEJAAAABuEJAAAABuEJAAAABuEJAAAABuEJAAAABvtCkmLFy9WdHS0/P4vP/isoqJCbrdbo0aNCtg2Ly9PLpdLxcXFVlt+fr5CQ0OVlZXVqvEWLFigoUOHKioqSr1799bYsWP15ptvtmcXAAAAbLUrJGVmZqqiokKbN2+22vLy8hQXF6fCwkIdPXrUavf5fIqPj1daWprVlpOTo7lz52r9+vXau3dvi+OlpaXp0Ucf1Y4dO7R+/XoNHjxY48aN0xdffNGe3QAAAGiiXSEpPT1d8fHx8vl8VpvP59OECROUkpKi/Pz8gPbMzExrubKyUsuXL9esWbN03XXX6cknn2xxvFtvvVVjx45VcnKyzj//fD388MMqKyvT9u3b27MbAAAATbT7NUkZGRnKzc21lnNzc5WRkSGv12u119bWauPGjQEh6dlnn1V6errS09M1depULVmyRG35M3K1tbV6/PHH5fF4NHz4cNttampqVFZWFvAAAABojQ4JSRs2bJDf71d5ebm2bt2q0aNHy+v1WneYCgoKVFVVFRCSsrOzNXXqVElSVlaWKioqtHbt2hbHW7VqlXr16qXIyEj97ne/02uvvaa+ffvabrtw4UJ5PB7rkZiY2N7dBQAA3US7Q1JmZqYqKytVWFiovLw8paWlKTY2Vl6vV4WFhaqsrJTP59OgQYOUnJwsSdq5c6c2bdqkKVOmSJLCwsI0efJk5eTktGq8oqIi5efnKysrS9/61rd08OBB223nz5+v0tJS67Fv37727i4AAOgmwtrbQWpqqhISEpSbm6uSkhJ5vV5JUlxcnIYMGaINGzYoNzdXY8aMsWqys7Pl9/s1cOBAq80YI7fbrZKSEvXu3TvoeFFRUUpNTVVqaqouu+wynXPOOcrOztb8+fObbBsREaGIiIj27iIAAOiGOuRzkjIzM+Xz+eTz+ZSRkWG1e71evfrqqyooKLB+1eb3+7V06VItWrRIRUVF1mPbtm1KSkrSsmXL2jS2MUY1NTUdsRsAAACWdt9Jko6HpNmzZ+vYsWPWnSTpeEiaNWuWqqurrZC0atUqlZSUaMaMGfJ4PAH93HTTTcrOztacOXOajFFZWalf/OIXuuGGGzRgwAAdOnRIf/7zn/XJJ5/o5ptv7ojdAAAAsHTYnaSqqiqlpqaqf//+VrvX61V5eblSUlKsF01nZ2dr7NixTQKSJE2aNElFRUXasmVLk3WhoaF6//33NWnSJKWlpem6667TF198oby8PJ1//vkdsRsAAACWDrmTNHjwYNu37yckJDRpX7lyZdB+LrnkkqAfAxAZGakVK1a0b6IAAACtxN9uAwAAsEFIAgAAsEFIAgAAsEFIAgAAsEFIAgAAsOEybfmrsme4srIyeTwelZaWKiYmpqun0yV8+59yVFdc6uyNkNV1jsoce2u/23HtuCHHHNWFO/yvRojLWd1ZEX5HdZu+cHZsDlc7KpMkffKFs0+8D49wduEMPMvZsfn4i3BHdQc+jXJUFxLi/MfuLVeWOKob1MvZMT1SE+qorszZt5P+tSvSWaGkw/9ydo3X9HD28+3sEc4+yHjqpZWO6r45ZKqjuq+Crnr+5k4SAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADZcxxnT1JDpLWVmZPB6PSktLFRMT09XT+cpYWPS/jur2lDsbLzbSWZ0k7TnkdlR3wYBaR3UxzoZTtMO6L6qd1YW4nNVJ0odlzorfefcsR3Xubccc1fndzv5PGHas3lFdeVoPR3WSdO2ofzuq6x3u7Mf53gpn5zAyzFGZYx/+2/mAFWXOvqkO/dvZefT9vyxHdbDXVc/f3EkCAACwQUgCAACwQUgCAACwQUgCAACwQUgCAACwQUgCAACwQUgCAACwQUgCAACw0a6QtHjxYkVHR8vv91ttFRUVcrvdGjVqVMC2eXl5crlcKi4uttry8/MVGhqqrKzWfejW9OnT5XK5Ah6XXXZZe3YBAADAVrtCUmZmpioqKrR582arLS8vT3FxcSosLNTRo0etdp/Pp/j4eKWlpVltOTk5mjt3rtavX6+9e/e2asysrCwdOHDAeqxevbo9uwAAAGCrXSEpPT1d8fHx8vl8VpvP59OECROUkpKi/Pz8gPbMzExrubKyUsuXL9esWbN03XXX6cknn2zVmBEREYqLi7Meffr0CbptTU2NysrKAh4AAACt0e7XJGVkZCg3N9dazs3NVUZGhrxer9VeW1urjRs3BoSkZ599Vunp6UpPT9fUqVO1ZMkStebPyPl8PsXGxiotLU0zZ87UwYMHg267cOFCeTwe65GYmNiOPQUAAN1Jh4SkDRs2yO/3q7y8XFu3btXo0aPl9XqtO0wFBQWqqqoKCEnZ2dmaOnWqpOO/QquoqNDatWubHWv8+PFatmyZXn/9dS1atEiFhYUaM2aMampqbLefP3++SktLrce+ffvau7sAAKCbaPffcM7MzFRlZaUKCwtVUlKitLQ0xcbGyuv16vbbb1dlZaV8Pp8GDRqk5ORkSdLOnTu1adMmrVix4vgkwsI0efJk5eTkaOzYsUHHmjx5svXvYcOGacSIEUpKStLLL7+siRMnNtk+IiJCERER7d1FAADQDbU7JKWmpiohIUG5ubkqKSmR1+uVJMXFxWnIkCHasGGDcnNzNWbMGKsmOztbfr9fAwcOtNqMMXK73SopKVHv3r1bNfaAAQOUlJSkXbt2tXc3AAAAAnTI5yRlZmbK5/PJ5/MpIyPDavd6vXr11VdVUFBg/arN7/dr6dKlWrRokYqKiqzHtm3blJSUpGXLlrV63EOHDmnfvn0aMGBAR+wGAACApd13kqTjIWn27Nk6duyYdSdJOh6SZs2aperqaiskrVq1SiUlJZoxY4Y8Hk9APzfddJOys7M1Z86cJmNUVFRowYIFmjRpkgYMGKA9e/bovvvuU9++ffXNb36zI3YDAADA0mF3kqqqqpSamqr+/ftb7V6vV+Xl5UpJSbHeWZadna2xY8c2CUiSNGnSJBUVFWnLli1N1oWGhmrHjh2aMGGC0tLSNG3aNKWlpWnjxo2Kjo7uiN0AAACwdMidpMGDB9u+fT8hIaFJ+8qVK4P2c8kllwT9GIAePXro1Vdfbd9EAQAAWom/3QYAAGCDkAQAAGCDkAQAAGCDkAQAAGCDkAQAAGDDZVrzV2W/IsrKyuTxeFRaWqqYmJiung6acd+m/3VUFx3ufMyPSl2O6gbHOPsWGhTlrK7O4XdsqLPd09slDgsl7S9xO6orLjzLUd3Qtz5zVOfUFwOdffxI1Nh6x2P2jPI7qot2WHd2pKMyhTv8L/jAns6fkt474uxa/dXXb3U8JjpHVz1/cycJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADABiEJAADAhssYY7p6Ep2lrKxMHo9HpaWliomJ6erpoAss27XMUd1nVc7GOzvCWV11nbO6EJezuiO1zuokaU+5s0GL3z/LUV31x87+bzfknX87qqs4y9lJ/DS1t6M6SYqJdXZCPGfVOKpLGujsAh8c7ezpIzLUUZnuPPc2Z4U443XV8zd3kgAAAGwQkgAAAGwQkgAAAGwQkgAAAGwQkgAAAGwQkgAAAGwQkgAAAGwQkgAAAGy0KyQtXrxY0dHR8vv9VltFRYXcbrdGjRoVsG1eXp5cLpeKi4uttvz8fIWGhiorK6tV47lcLtvHb37zm/bsBgAAQBPtCkmZmZmqqKjQ5s2brba8vDzFxcWpsLBQR48etdp9Pp/i4+OVlpZmteXk5Gju3Llav3699u7d2+J4Bw4cCHjk5OTI5XJp0qRJ7dkNAACAJtoVktLT0xUfHy+fz2e1+Xw+TZgwQSkpKcrPzw9oz8zMtJYrKyu1fPlyzZo1S9ddd52efPLJFseLi4sLeLz44ovKzMxUcnKy7fY1NTUqKysLeAAAALRGu1+TlJGRodzcXGs5NzdXGRkZ8nq9Vnttba02btwYEJKeffZZpaenKz09XVOnTtWSJUvUlj8j9/nnn+vll1/WjBkzgm6zcOFCeTwe65GYmOhgDwEAQHfUISFpw4YN8vv9Ki8v19atWzV69Gh5vV7rDlNBQYGqqqoCQlJ2dramTp0qScrKylJFRYXWrl3b6nGfeuopRUdHa+LEiUG3mT9/vkpLS63Hvn37nO0kAADodtodkjIzM1VZWanCwkLl5eUpLS1NsbGx8nq9KiwsVGVlpXw+nwYNGmT9Wmznzp3atGmTpkyZIkkKCwvT5MmTlZOT0+pxc3JydNtttykyMjLoNhEREYqJiQl4AAAAtEZYeztITU1VQkKCcnNzVVJSIq/XK+n464eGDBmiDRs2KDc3V2PGjLFqsrOz5ff7NXDgQKvNGCO3262SkhL17t272THz8vK0c+dOPfvss+2dPgAAgK0O+ZykzMxM+Xw++Xw+ZWRkWO1er1evvvqqCgoKrF+1+f1+LV26VIsWLVJRUZH12LZtm5KSkrRs2bIWx8vOztall16q4cOHd8T0AQAAmmj3nSTpeEiaPXu2jh07Zt1Jko6HpFmzZqm6utoKSatWrVJJSYlmzJghj8cT0M9NN92k7OxszZkzJ+hYZWVleu6557Ro0aKOmDoAAICtDruTVFVVpdTUVPXv399q93q9Ki8vV0pKivXOsuzsbI0dO7ZJQJKkSZMmqaioSFu2bAk61jPPPCNjjG655ZaOmDoAAICtDrmTNHjwYNu37yckJDRpX7lyZdB+LrnkkhY/BuDOO+/UnXfe6WyiAAAArcTfbgMAALBBSAIAALBBSAIAALBBSAIAALBBSAIAALDhMm35q7JnuLKyMnk8HpWWlvInStBhnv+o5Q9AtfPpUWfjHal1Oaqrc/idXu13VidJ/6509v+w3R80/YiQ1ojYUeuozlXv7OAcTIx2VBcW4fzH7sDEckd1fftXOR7Tib+MurVTx8NXW1c9f3MnCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwIbLGGO6ehKdpaysTB6PR6WlpYqJienq6QBB/eHtZY7qDtW4HNWFOiuTJB2sclZXciTcUd3uDzyO6lyH6x3V1US5HdUNOeeIozpJ6hnld1QXHlHnqO4J7y2O6oDO0lXP39xJAgAAsEFIAgAAsEFIAgAAsEFIAgAAsEFIAgAAsEFIAgAAsEFIAgAAsEFIAgAAsNGukLR48WJFR0fL7//yg88qKirkdrs1atSogG3z8vLkcrlUXFxsteXn5ys0NFRZWVmtGu/zzz/X9OnTFR8fr549eyorK0u7du1qzy4AAADYaldIyszMVEVFhTZv3my15eXlKS4uToWFhTp69KjV7vP5FB8fr7S0NKstJydHc+fO1fr167V3795mxzLG6MYbb9RHH32kF198UVu3blVSUpLGjh2rysrK9uwGAABAE+0KSenp6YqPj5fP57PafD6fJkyYoJSUFOXn5we0Z2ZmWsuVlZVavny5Zs2apeuuu05PPvlks2Pt2rVLBQUFeuyxx/S1r31N6enp+vOf/6yKigo9/fTTtjU1NTUqKysLeAAAALRGu1+TlJGRodzcXGs5NzdXGRkZ8nq9Vnttba02btwYEJKeffZZpaenKz09XVOnTtWSJUvU3J+Rq6mpkSRFRkZabaGhoQoPD9f69ettaxYuXCiPx2M9EhMT27WvAACg++iQkLRhwwb5/X6Vl5dr69atGj16tLxer3WHqaCgQFVVVQEhKTs7W1OnTpUkZWVlqaKiQmvXrg06ztChQ5WUlKT58+erpKREtbW1euihh/TZZ5/pwIEDtjXz589XaWmp9di3b197dxcAAHQT7Q5JmZmZqqysVGFhofLy8pSWlqbY2Fh5vV4VFhaqsrJSPp9PgwYNUnJysiRp586d2rRpk6ZMmSJJCgsL0+TJk5WTkxN0HLfbreeff17FxcXq06ePevbsKZ/Pp/Hjxys0NNS2JiIiQjExMQEPAACA1ghrbwepqalKSEhQbm6uSkpK5PV6JUlxcXEaMmSINmzYoNzcXI0ZM8aqyc7Olt/v18CBA602Y4zcbrdKSkrUu3dv27EuvfRSFRUVqbS0VLW1terXr5++/vWva8SIEe3dDQAAgAAd8jlJmZmZ8vl88vl8ysjIsNq9Xq9effVVFRQUWL9q8/v9Wrp0qRYtWqSioiLrsW3bNiUlJWnZsmUtjufxeNSvXz/t2rVLmzdv1oQJEzpiNwAAACztvpMkHQ9Js2fP1rFjx6w7SdLxkDRr1ixVV1dbIWnVqlUqKSnRjBkz5PF4Avq56aablJ2drTlz5tiO89xzz6lfv34aNGiQduzYoXnz5unGG2/UuHHjOmI3AAAALB12J6mqqkqpqanq37+/1e71elVeXq6UlBTrnWXZ2dkaO3Zsk4AkSZMmTVJRUZG2bNliO86BAwd0++23a+jQobrrrrt0++23B337PwAAQHu4THPvu/+KKSsrk8fjUWlpKS/ixmntD2+3/GtnO4dqXI7qQp2VSZIOVjmrKzkS7qhu9wdN/4PVGq7D9Y7qaqLcjuqGnHPEUZ0k9Yzyt7yRjfCIOkd1T3hvcVQHdJauev7mb7cBAADYICQBAADYICQBAADYICQBAADYICQBAADY4N1tQDfh9B1zR2qdv/XtULWzuk/293BUV3okwlHd2f2cTbRXdK2jujC3s3faSdJfRt3quBY4U/HuNgAAgNMIIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMAGIQkAAMCGyxhjunoSnaWsrEwej0elpaWKiYnp6ukAXwkLi/7XUd3hamfjVdU5q4sKc1ZX7/An5G8uu9VZIYAmuur5mztJAAAANghJAAAANghJAAAANghJAAAANghJAAAANghJAAAANghJAAAANghJAAAANlodkhYvXqzo6Gj5/X6rraKiQm63W6NGjQrYNi8vTy6XS8XFxZKk/Px8hYaGKisry7bvefPm6dJLL1VERIQuuuiiJuurq6s1ffp0XXDBBQoLC9ONN97Y2mkDAAA40uqQlJmZqYqKCm3evNlqy8vLU1xcnAoLC3X06FGr3efzKT4+XmlpaZKknJwczZ07V+vXr9fevXub9G2M0R133KHJkyfbjl1XV6cePXrorrvu0tixY1u9cwAAAE61OiSlp6crPj5ePp/PavP5fJowYYJSUlKUn58f0J6ZmSlJqqys1PLlyzVr1ixdd911evLJJ5v0/Yc//EGzZ89WcnKy7dhRUVF67LHHNHPmTMXFxbV2yqqpqVFZWVnAAwAAoDXa9JqkjIwM5ebmWsu5ubnKyMiQ1+u12mtra7Vx40YrJD377LNKT09Xenq6pk6dqiVLlqiz/lzcwoUL5fF4rEdiYmKnjAsAAM58bQ5JGzZskN/vV3l5ubZu3arRo0fL6/Vad5gKCgpUVVVlhaTs7GxNnTpVkpSVlaWKigqtXbu2Y/ciiPnz56u0tNR67Nu3r1PGBQAAZ742haTMzExVVlaqsLBQeXl5SktLU2xsrLxerwoLC1VZWSmfz6dBgwYpOTlZO3fu1KZNmzRlyhRJUlhYmCZPnqycnJxTsjONRUREKCYmJuABAADQGmFt2Tg1NVUJCQnKzc1VSUmJvF6vJCkuLk5DhgzRhg0blJubqzFjxkg6fhfJ7/dr4MCBVh/GGLndbpWUlKh3794duCsAAAAdp82fk5SZmSmfzyefz6eMjAyr3ev16tVXX1VBQYEyMzPl9/u1dOlSLVq0SEVFRdZj27ZtSkpK0rJlyzpyPwAAADpUm+4kScdD0uzZs3Xs2DHrTpJ0PCTNmjVL1dXVyszM1KpVq1RSUqIZM2bI4/EE9HHTTTcpOztbc+bMkSR98MEHqqio0GeffaaqqioVFRVJks477zyFh4dLkt59913V1tbq8OHDKi8vt7ax+1wlAACA9nKZNr7VbM+ePRoyZIiGDh2q9957z2r/5JNPlJiYqJSUFH3wwQe6/vrrVV9fr5dffrlJH1u2bNGll16qt956S5dccokyMjK0bt26Jtvt3r1bgwcPliQNHjxYH3/8cZNt2jL9srIyeTwelZaW8vokoIMsLPpfR3WHq52NV1XnrC6qzf8lPK7e4Ztxf3PZrc4KATTRVc/fbf6xMXjwYNtgkpCQENC+cuXKoH1ccsklAdue/NlLwezZs6dN8wQAAGgP/nYbAACADUISAACADUISAACADUISAACADUISAACADYdvigWA4+Zf1PFvdX/mA2cfNnuoxtl4s8+/zVkhgK807iQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYCOvqCXQmY4wkqaysrItnAqA5R8uPOqqrqnU2Hj8TgNNbw/dow/N4Z3GZzh6xC33yySdKTEzs6mkAAAAH9u3bp4SEhE4br1uFpPr6eu3fv1/R0dFyuVxdPZ1ToqysTImJidq3b59iYmK6ejqnFY6NPY5LcByb4Dg2wXFsgnN6bIwxKi8vV3x8vEJCOu+VQt3q120hISGdmkC7UkxMDN+cQXBs7HFcguPYBMexCY5jE5yTY+PxeE7RbILjhdsAAAA2CEkAAAA2CElfMREREXrwwQcVERHR1VM57XBs7HFcguPYBMexCY5jE9yZdmy61Qu3AQAAWos7SQAAADYISQAAADYISQAAADYISQAAADYISV8BCxYskMvlCnjExcV19bS6xBtvvKHrr79e8fHxcrlc+uc//xmw3hijBQsWKD4+Xj169FBGRobeeeedrplsJ2vp2EyfPr3JdXTZZZd1zWQ70cKFC/W1r31N0dHRio2N1Y033qidO3cGbNNdr5vWHJvuet089thjuvDCC60PRRw5cqReeeUVa313vWaklo/NmXTNEJK+Is4//3wdOHDAeuzYsaOrp9QlKisrNXz4cD366KO263/961/r4Ycf1qOPPqrCwkLFxcXp6quvVnl5eSfPtPO1dGwkKSsrK+A6Wr16dSfOsGusW7dOs2fPVkFBgV577TX5/X6NGzdOlZWV1jbd9bppzbGRuud1k5CQoIceekibN2/W5s2bNWbMGE2YMMEKQt31mpFaPjbSGXTNGJzxHnzwQTN8+PCunsZpR5J54YUXrOX6+noTFxdnHnroIauturraeDwes3jx4i6YYddpfGyMMWbatGlmwoQJXTKf08nBgweNJLNu3TpjDNfNyRofG2O4bk7Wu3dv89e//pVrxkbDsTHmzLpmuJP0FbFr1y7Fx8dryJAhmjJlij766KOuntJpZ/fu3frss880btw4qy0iIkJer1f5+fldOLPTh8/nU2xsrNLS0jRz5kwdPHiwq6fU6UpLSyVJffr0kcR1c7LGx6ZBd79u6urq9Mwzz6iyslIjR47kmjlJ42PT4Ey5ZrrVH7j9qvr617+upUuXKi0tTZ9//rl+/vOf6/LLL9c777yjs88+u6und9r47LPPJEn9+/cPaO/fv78+/vjjrpjSaWX8+PG6+eablZSUpN27d+u///u/NWbMGL311ltnzKfjtpcxRt///vd15ZVXatiwYZK4bhrYHRupe183O3bs0MiRI1VdXa1evXrphRde0HnnnWcFoe58zQQ7NtKZdc0Qkr4Cxo8fb/37ggsu0MiRI5WSkqKnnnpK3//+97twZqcnl8sVsGyMadLWHU2ePNn697BhwzRixAglJSXp5Zdf1sSJE7twZp1nzpw52r59u9avX99kXXe/boIdm+583aSnp6uoqEhHjhzR888/r2nTpmndunXW+u58zQQ7Nuedd94Zdc3w67avoKioKF1wwQXatWtXV0/ltNLwjr+GOwMNDh482OR/fJAGDBigpKSkbnMdzZ07Vy+99JJyc3OVkJBgtXPdBD82drrTdRMeHq7U1FSNGDFCCxcu1PDhw/X73/+ea0bBj42d0/maISR9BdXU1Oi9997TgAEDunoqp5UhQ4YoLi5Or732mtVWW1urdevW6fLLL+/CmZ2eDh06pH379n3lryNjjObMmaMVK1bo9ddf15AhQwLWd+frpqVjY6e7XDd2jDGqqanp1tdMMA3Hxs5pfc102UvG0WF+8IMfGJ/PZz766CNTUFBgrrvuOhMdHW327NnT1VPrdOXl5Wbr1q1m69atRpJ5+OGHzdatW83HH39sjDHmoYceMh6Px6xYscLs2LHD3HLLLWbAgAGmrKysi2d+6jV3bMrLy80PfvADk5+fb3bv3m1yc3PNyJEjzcCBA7/yx2bWrFnG4/EYn89nDhw4YD2OHj1qbdNdr5uWjk13vm7mz59v3njjDbN7926zfft2c99995mQkBCzZs0aY0z3vWaMaf7YnGnXDCHpK2Dy5MlmwIABxu12m/j4eDNx4kTzzjvvdPW0ukRubq6R1OQxbdo0Y8zxt3M/+OCDJi4uzkRERJjRo0ebHTt2dO2kO0lzx+bo0aNm3Lhxpl+/fsbtdptBgwaZadOmmb1793b1tE85u2MiySxZssTaprteNy0dm+583dxxxx0mKSnJhIeHm379+pmrrrrKCkjGdN9rxpjmj82Zds24jDGm8+5bAQAAnBl4TRIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAIANQhIAAICNsK6eQGerrq5WbW1tV08DAAC0QXh4uCIjIzt1zG4Vkqqrq9XnrDhV1ZR29VQAAEAbxMXFaffu3Z0alLpVSKqtrVVVTakmjXtEoeE9jzee+IVjfajr+NeQwK/G5bLWf7nueJEJcdnWqqHdZb/ehNiPVR+qwPWuL7c3QcYICTGSJFeoCVgOcZ342rj9xP5ay43WuxrVh4aak2ob1zTqq3Efwb42GjM0cJesr6GuputOHKIv20Oa1gQuB47R3FjBxm746go2h6B9mUbr7Zeb9mOa9OFqtD/WtjqxrCDrW1huvL1LIXKdGMx1otcmy2rYuGG54fgErm+pvsn2rhCbbVzN1zQa48ThkEx9kK8nNqi3Wd/wb7t1J301jbdT89sH7++kuQSrqTfNtzdZDtbeeB9Oam9zn0GOZZB+TF2jfk4sG2OkOp1oO7Fto20aaqw+Go5DQx+Nt2/S3mh/T+rXNIzdeNtG+/3lqXE12n2X7fqgy/pyuekl0Whb02isxtuZ5uuCzbm+zlj/rjvRSf2J41BvHevA5ZNrTzpMTdc3jNmovu7ksc2X/w6sPbFtXcO8m86tSn7d81m+amtrCUmnmjush8LcPSQFDzpNAk0bQlKTENTqkBR8fUNIcgUJSU3DUAvtrVzvJCS1dUy7gNLwtbl1J88hWAhqKSQ1197Wmqahp/lA0paQ1FIfLYWk1vbXlpBkF3ICljskJLVtjJZDUjPBpdUhqa5V27W2v1aFpBaCSMtBJ3C71oWkjmkPGpLqzUnhppUh6ctnePvtg4Ukm36bzst+27aGpPpGy3YBpr6lQNVSSLIOeRtDkuukkNQQRPTlOtvlE99a9Se+ueoaLX/5X4UTy6bR9g3rjQmYR0DfrkZ9261v+N7uZLxwGwAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwAYhCQAAwEZYV0+gKxzzV6k+xHV84URMrK8/vtzQ3vDVuE4sh7pOWne8yIR8ue7kGjW0u+zXmxD7sepDFbje9eX2JsgYJsScqD3+NeTEcojrxNfG7Sf211putN7VqL4+1JxU27imUV+N+wj2tdGYoYG7ZH0NdTVdd+IQfdke0rQmcDlwjObGCjZ2w1dXsDkE7cs0Wm+/3LQf06QPV6P9+fLyPbGsIOtbWG68vUshcp0YzHWi1ybLati4Ybnh+ASub6m+yfauEJttXM3XNBrjxOGQTH2Qryc2qLdZ3/Bvu3UnfTWNt1Pz2wfv76S5BKupN823N1kO1t54H05qb3OfQY5lkH5MXaN+TiwbY6Q6nWg7sW2jbRpqrD4ajkNDH423b9LeaH9P6tc0jN1420b7/eWpcTXafZft+qDL+nK56SXRaFvTaKzG25nm64LNub7OWP+uO9FJ/YnjUG8d68Dlk2tPOkxN1ze5hJtu/+WhPVHbaNk6zTbrq+RXV3AZ6+r56istLVVCQoIqKiq6eioAAKANevXqpU8++UQej6fTxuxWd5JcLpcqKiq0b98+xcTEdPV0cIqVlZUpMTGR891NcL67F85399Jwvq27xp2kW4WkBjExMXxTdSOc7+6F8929cL5xKvHCbQAAABuEJAAAABvdKiRFRETowQcfVERERFdPBZ2A8929cL67F85399JV57tbvbsNAACgtbrVnSQAAIDWIiQBAADYICQBAADYICQBAADYOGND0qeffqqpU6fq7LPPVs+ePXXRRRfprbfearZm3bp1uvTSSxUZGank5GQtXrw4YP2KFSs0YsQInXXWWYqKitJFF12kv/3tb6dyNxDEG2+8oeuvv17x8fFyuVz65z//GbDeGKMFCxYoPj5ePXr0UEZGht55550W+33kkUeUnp6uHj16KDExUXfffbeqq6ut9X6/Xw888ICGDBmiHj16KDk5WT/96U9Vb/2dLpwKLZ3vzz//XNOnT1d8fLx69uyprKws7dq1q8V+jxw5otmzZ2vAgAGKjIzUueeeq9WrV1vrH3vsMV144YXWBxKOHDlSr7zySkfvHlrQ0nmYPn26XC5XwOOyyy5rdf/PPPOMXC6XbrzxxlMwe7SkpefrBQsWaOjQoYqKilLv3r01duxYvfnmmy3229LP85MtXLhQLpdL//Vf/9WmuZ+Rn7hdUlKiK664QpmZmXrllVcUGxurDz/8UGeddVbQmt27d+sb3/iGZs6cqb///e/asGGDvve976lfv36aNGmSJKlPnz66//77NXToUIWHh2vVqlX6zne+o9jYWF1zzTWdtHeQpMrKSg0fPlzf+c53rPNzsl//+td6+OGH9eSTTyotLU0///nPdfXVV2vnzp2Kjo627XPZsmW69957lZOTo8svv1zFxcWaPn26JOl3v/udJOlXv/qVFi9erKeeekrnn3++Nm/erO985zvyeDyaN2/eKdvf7q65822M0Y033ii3260XX3xRMTExevjhhzV27Fi9++67ioqKsu2ztrZWV199tWJjY/WPf/xDCQkJ2rdvX8D1kZCQoIceekipqamSpKeeekoTJkzQ1q1bdf7555+6HUaA1pyHrKwsLVmyxKoJDw9vVd8ff/yx7rnnHo0aNarjJ44Wteb5Oi0tTY8++qiSk5NVVVWl3/3udxo3bpw++OAD9evXz7bf1vw8b1BYWKjHH39cF154Ydt3wJyBfvzjH5srr7yyTTU/+tGPzNChQwPa/vM//9NcdtllzdZdfPHF5oEHHmjzHNFxJJkXXnjBWq6vrzdxcXHmoYcestqqq6uNx+MxixcvDtrP7NmzzZgxYwLavv/97wdcS9dee6254447AraZOHGimTp1ajv3Aq3V+Hzv3LnTSDJvv/221eb3+02fPn3ME088EbSfxx57zCQnJ5va2to2jd+7d2/z17/+tc3zRsc6+TxMmzbNTJgwoc19+P1+c8UVV5i//vWvjvtA+zh5vi4tLTWSzL/+9a+g27Tm57kxxpSXl5tzzjnHvPbaa8br9Zp58+a1aS5n5K/bXnrpJY0YMUI333yzYmNjdfHFF+uJJ55otmbjxo0aN25cQNs111yjzZs369ixY022N8Zo7dq12rlzp0aPHt2h80f77N69W5999lnA+YyIiJDX61V+fn7QuiuvvFJvvfWWNm3aJEn66KOPtHr1al177bUB26xdu1bFxcWSpG3btmn9+vX6xje+cYr2Bi2pqamRJEVGRlptoaGhCg8P1/r164PWvfTSSxo5cqRmz56t/v37a9iwYfrlL3+puro62+3r6ur0zDPPqLKyUiNHjuzYnUCrBTsPPp9PsbGxSktL08yZM3Xw4MEW+/rpT3+qfv36acaMGadyymhGW5+va2tr9fjjj8vj8Wj48OFBt2vNz3NJmj17tq699lqNHTvW2Q60KVKdJiIiIkxERISZP3++2bJli1m8eLGJjIw0Tz31VNCac845x/ziF78IaNuwYYORZPbv32+1HTlyxERFRZmwsDATERFhsrOzT9l+oHXU6M5Cw3n79NNPA7abOXOmGTduXLN9/eEPfzBut9uEhYUZSWbWrFkB6+vr6829995rXC6XCQsLMy6Xy/zyl7/ssH1Byxqf79raWpOUlGRuvvlmc/jwYVNTU2MWLlxoJDV7vtPT001ERIS54447zObNm83TTz9t+vTpY37yk58EbLd9+3YTFRVlQkNDjcfjMS+//PKp2jU0o7nz8Mwzz5hVq1aZHTt2mJdeeskMHz7cnH/++aa6ujpof+vXrzcDBw40X3zxhTHG+d0otE9rn69XrlxpoqKijMvlMvHx8WbTpk0t9t3Sz/Onn37aDBs2zFRVVRljjKM7SWdkSHK73WbkyJEBbXPnzm32V2fnnHNOkye79evXG0nmwIEDVltdXZ3ZtWuX2bp1q/ntb39rPB6Pyc3N7dD5o22ChaSTw60xxnz3u98111xzTdB+cnNzTf/+/c0TTzxhtm/fblasWGESExPNT3/6U2ubp59+2iQkJJinn37abN++3SxdutT06dPHPPnkkx2+X7DX+HwbY8zmzZvN8OHDjSQTGhpqrrnmGjN+/Hgzfvz4oP2cc845JjEx0fj9fqtt0aJFJi4uLmC7mpoas2vXLlNYWGjuvfde07dvX/POO+906D6hZW05D/v37zdut9s8//zztuvLysrM4MGDzerVq602QlLXaO3zdUVFhdm1a5fZuHGjueOOO8zgwYPN559/HrTfln6e792718TGxpqioiKrptuEpEGDBpkZM2YEtP35z3828fHxQWtGjRpl7rrrroC2FStWmLCwsGZfszBjxowW707g1Gr8pPnhhx8aSWbLli0B291www3m29/+dtB+rrzySnPPPfcEtP3tb38zPXr0MHV1dcYYYxISEsyjjz4asM3PfvYzk56e3s69QGvZhaQGR44cMQcPHjTGGPMf//Ef5nvf+17QfkaPHm2uuuqqgLbVq1cbSaampiZo3VVXXWXuvPPOtk8cHaql85CamhrwusSTbd261QrUDQ+Xy2VcLpcJDQ01H3zwwamaNhpx8nxtzPHz29xd/JZ+nr/wwgtNrgFJ1jVw8n+emnNGvibpiiuu0M6dOwPaiouLlZSUFLRm5MiReu211wLa1qxZoxEjRsjtdgetM8ZYr4nA6WHIkCGKi4sLOJ+1tbVat26dLr/88qB1R48eVUhI4CUfGhoqc/w/C81uw0cAnB48Ho/69eunXbt2afPmzZowYULQba+44gp98MEHAeeuuLhYAwYMaPadUXzPnx6aOw+HDh3Svn37NGDAANv1Q4cO1Y4dO1RUVGQ9brjhBmVmZqqoqEiJiYmncuo4iZPna6nl78OWfp5fddVVTa6BESNG6LbbblNRUZFCQ0NbtwOtilKnmU2bNpmwsDDzi1/8wuzatcssW7bM9OzZ0/z973+3trn33nvN7bffbi1/9NFHpmfPnubuu+827777rsnOzjZut9v84x//sLb55S9/adasWWM+/PBD895775lFixaZsLCwZt9Bg1OjvLzcbN261fof4cMPP2y2bt1qPv74Y2OMMQ899JDxeDxmxYoVZseOHeaWW24xAwYMMGVlZVYft99+u7n33nut5QcffNBER0ebp59+2nz00UdmzZo1JiUlxXzrW9+ytpk2bZoZOHCgWbVqldm9e7dZsWKF6du3r/nRj37UeTvfDbV0vpcvX25yc3PNhx9+aP75z3+apKQkM3HixIA+Gp/vvXv3ml69epk5c+aYnTt3mlWrVpnY2Fjz85//3Npm/vz55o033jC7d+8227dvN/fdd58JCQkxa9as6ZwdhzGm+fNQXl5ufvCDH5j8/Hyze/duk5uba0aOHGkGDhzY7Pd7Y/y6rWu09HxdUVFh5s+fbzZu3Gj27Nlj3nrrLTNjxgwTERER8I5WJz/PG+s2v24z5viLvIYNG2YiIiLM0KFDzeOPPx6wftq0acbr9Qa0+Xw+c/HFF5vw8HAzePBg89hjjwWsv//++01qaqqJjIw0vXv3NiNHjjTPPPPMqd4V2MjNzTWSmjymTZtmjDn+AusHH3zQxMXFmYiICDN69GizY8eOgD68Xq+1vTHGHDt2zCxYsMCkpKSYyMhIk5iYaL73ve+ZkpISa5uysjIzb948M2jQIBMZGWmSk5PN/fff3+yvZ9B+LZ3v3//+9yYhIcG43W4zaNAg88ADDzQ5J43PtzHG5Ofnm69//esmIiLCJCcnm1/84hcBt9nvuOMOk5SUZMLDw02/fv3MVVddRUDqAs2dh6NHj5px48aZfv36Wed/2rRpZu/evQF92J3/kxGSuk5zz9dVVVXmm9/8pomPjzfh4eFmwIAB5oYbbmjywm0nP88bcxKSXMac+D0DAAAALGfka5IAAABONUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISAACADUISgK+sjIwM/dd//Ven1bd3vOYcOnRIsbGx2rNnzynp385NN92khx9+uNPGA043/O02oAvl5+dr1KhRuvrqq/V///d/XT2dDpeRkaGLLrpIjzzySJeMf/jwYbndbkVHR7e4rd1c21Jvt31H7v8999yjkpISZWdnB7Rv375dCxYs0Lp161RZWamUlBRNnz5dd999t8LCwto15vbt25WZmandu3crJiamXX0BZyLuJAFdKCcnR3PnztX69eu1d+/erp6OJKm2trarp9Bhc+jTp0+rA05H1Ld3vGCqqqqUnZ2t7373uwHt69at02WXXaYePXroxRdf1LZt2/SjH/1Iv/3tbzVx4kTV19e3a9wLL7xQgwcP1rJly9rVD3DGMgC6REVFhYmOjjbvv/++mTx5svnJT34SsN7r9ZrZs2eb2bNnG4/HY/r06WPuv/9+U19f36ZtXnnlFXPFFVdY66+99lrzwQcfNOnj7rvvNmeffbYZPXq0McaY+vp686tf/coMGTLEREZGmgsvvNA899xzAXVz5841P/zhD03v3r1N//79zYMPPmitnzZtmpEU8Ni9e7ftsXA6h+eee84MGzbMREZGmj59+pirrrrKVFRUBPQ7b968FrcNNteT6xcvXmzi4+NNXV1dwNyvv/568+1vf7vJeHZ9/uQnPzF9+vQx1dXVAX1MnDjR3H777bbHxhhjnn/+edO3b9+ANr/fb4YMGWJuu+22Jtu/8847xu12m7/+9a/WvObMmWPmzZtnzjrrLBMbG2v+8pe/mIqKCjN9+nTTq1cvk5ycbFavXt2krwULFphRo0YFnRvwVUZIArpIdna2GTFihDHGmJUrV5rBgwc3CUC9evUy8+bNM++//775+9//bnr27Gkef/zxNm3zj3/8wzz//POmuLjYbN261Vx//fXmggsusJ7sG/r44Q9/aN5//33z3nvvGWOMue+++8zQoUPN//3f/5kPP/zQLFmyxERERBifz2fVxcTEmAULFpji4mLz1FNPGZfLZdasWWOMMebIkSNm5MiRZubMmebAgQPmwIEDxu/32x4LJ3PYv3+/CQsLMw8//LDZvXu32b59u/nTn/5kysvLA/qdN29ei9sGm+vJoefQoUMmPDzc/Otf/7L6P3z4sAkPDzevvvpqwHjB+qyoqDAej8csX77c6uOLL74w4eHh5vXXXw96rcybN89kZWUFtOXn5xtJpqioyLZmwoQJ5uqrr7bmFR0dbX72s5+Z4uJi87Of/cyEhISY8ePHm8cff9wUFxebWbNmmbPPPttUVlYG9LN69WoTERHRJNgB3QEhCegil19+uXnkkUeMMcYcO3bM9O3b17z22mvWeq/Xa84999yA4PTjH//YnHvuuW3aprGDBw8aSWbHjh1WHxdddFHANhUVFSYyMtLk5+cHtM+YMcPccsstVt2VV14ZsP5rX/ua+fGPfxwwv4bQ0Bwnc3jrrbeMJLNnz55m+503b16btm2u7YYbbjB33HGHtfyXv/zFxMXFWeGv8fZ2fc6aNcuMHz/eWn7kkUdMcnJywDlsbMKECQHjGmPMsmXLjKQmoabBD3/4Q5OWlmbN4+Rz5ff7TVRUVMDdqwMHDhhJZuPGjQH9bNu2rcVjB3xV8ZokoAvs3LlTmzZt0pQpUyRJYWFhmjx5snJycgK2u+yyy+RyuazlkSNHateuXaqrq2v1Nh9++KFuvfVWJScnKyYmRkOGDJGkgNdAjRgxImDcd999V9XV1br66qvVq1cv67F06VJ9+OGH1nYXXnhhQN2AAQN08OBBR8ekrXMYPny4rrrqKl1wwQW6+eab9cQTT6ikpMS277Zs25zbbrtNzz//vGpqaiRJy5Yt05QpUxQaGtrqPmbOnKk1a9bo008/lSQtWbJE06dPDziHjVVVVSkyMjKgreGF1IcPH7atKSkpCXix9cnnKjQ0VGeffbYuuOACq61///6S1OT89ejRQ5J09OjRFvcN+Kpp31sfADiSnZ0tv9+vgQMHWm3GGLndbpWUlKh3794dNtb111+vxMREPfHEE4qPj1d9fb2GDRsW8OLoqKiogJqGF/y+/PLLAXOUpIiICOvfbrc7YJ3L5XL8YuG2ziE0NFSvvfaa8vPztWbNGv3xj3/U/fffrzfffNMKgg3asm1zrr/+etXX1+vll1/W1772NeXl5bX5LfIXX3yxhg8frqVLl+qaa67Rjh07tHLlymZr+vbt2yTUjRw5Um63WytXrtSsWbMC1tXV1WnNmjWaOHGi1WZ3rk5uawhpjc9fQwjr169fK/cQ+OrgThLQyfx+v5YuXapFixapqKjIemzbtk1JSUkB7yQqKCgIqC0oKNA555wTcOeiuW0OHTqk9957Tw888ICuuuoqnXvuua26g3LeeecpIiJCe/fuVWpqasAjMTGx1fsaHh4ecNerLVozB5fLpSuuuEI/+clPtHXrVoWHh+uFF16w7a+lbVsz1x49emjixIlatmyZnn76aaWlpenSSy8Nun2wPr/73e9qyZIlysnJ0dixY1s8phdffLHefffdgLazzz5bd911l37+859r//79Aet+97vf6dChQ7r77rub7bc13n77bSUkJKhv377t7gs403AnCehkq1atUklJiWbMmCGPxxOw7qabblJ2drbmzJkjSdq3b5++//3v6z//8z+1ZcsW/fGPf9SiRYsCaprbpnfv3jr77LP1+OOPa8CAAdq7d6/uvffeFucYHR2te+65R3fffbfq6+t15ZVXqqysTPn5+erVq5emTZvWqn0dPHiw3nzzTe3Zs0e9evVSnz59FBLSuv+btTSHoUOHau3atRo3bpxiY2P15ptv6osvvtC5557bpK8333yzxW3t5mrntttu0/XXX6933nlHU6dOdbT/t912m+655x498cQTWrp0aYvH4pprrtH8+fMD7jJWVFTorrvuUkFBgTIzM/X000/rkksu0W9+8xvdf//9+stf/tKukNogLy9P48aNa1cfwJmKO0lAJ8vOztbYsWObBCRJmjRpkoqKirRlyxZJ0re//W1VVVXpP/7jPzR79mzNnTtXd955Z0BNc9uEhITomWee0VtvvaVhw4bp7rvv1m9+85tWzfNnP/uZ/ud//kcLFy7Uueeeq2uuuUYrV65s06+n7rnnHoWGhuq8885Tv3792vxZUM3NISYmRm+88Ya+8Y1vKC0tTQ888IAWLVqk8ePHN+mnNdu2dq5jxoxRnz59tHPnTt16662O9j8mJkaTJk1Sr169dOONN7Z4HC644AKNGDFCy5cvt9p++9vfKikpSRs2bFBxcbH+8Ic/SJJ+9KMf6dixY7rjjjs0YMAA7du3r8X+g6murtYLL7ygmTNnOu4DOJPxidvAaao1n9bc1Z9oDeeuvvpqnXvuuVa4acnq1at1zz336O2332713bj2+tOf/qQXX3xRa9as6ZTxgNMNv24DgE50+PBhrVmzRq+//roeffTRVtd94xvf0K5du/Tpp5+26XVh7eF2u/XHP/6xU8YCTkeEJADoRJdccolKSkr0q1/9Sunp6W2qnTdv3imalb3Gv9oFuht+3QYAAGCDF24DAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADYICQBAADY+P+ZZA6DxCBq1gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data2d = data.copy()  # pg.DataContainerERT(data)\n",
    "xa = pg.x(data2d)[data2d[\"a\"]]\n",
    "xm = pg.x(data2d)[data2d[\"m\"]]\n",
    "data2d[\"valid\"] = (xa == 220.) * (xm == 220.)\n",
    "data2d.removeInvalid()\n",
    "data2d.removeUnusedSensors()\n",
    "data2d.swapXY()  # exchange x and y so that profile runs along x\n",
    "data2d.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "28/08/23 - 11:14:24 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Found 2 regions.\n",
      "28/08/23 - 11:14:24 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - (ERTModelling) Region with smallest marker (1) set to background.\n",
      "28/08/23 - 11:14:24 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Creating forward mesh from region infos.\n",
      "28/08/23 - 11:14:26 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Creating refined mesh (H2) to solve forward task.\n",
      "28/08/23 - 11:14:29 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Mesh for forward task: Mesh: Nodes: 38336 Cells: 196784 Boundaries: 208012\n",
      "28/08/23 - 11:14:31 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Use median(data values)=11.029958607999644\n",
      "28/08/23 - 11:14:31 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Created startmodel from forward operator: 20703, min/max=11.029959/11.029959\n",
      "28/08/23 - 11:14:31 - pyGIMLi - \u001b[0;32;49mINFO\u001b[0m - Starting inversion.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x0000020151106310>\n",
      "Data transformation: <pgcore._pygimli_.RTransLogLU object at 0x000002013F3791D0>\n",
      "Model transformation (cumulative):\n",
      "\t 0 <pgcore._pygimli_.RTransLogLU object at 0x000002014EAB69A0>\n",
      "min/max (data): 4.84/201\n",
      "min/max (error): 5%/6.52%\n",
      "min/max (start model): 11.03/11.03\n",
      "--------------------------------------------------------------------------------\n",
      "--------------------------------------------------------------------------------\n",
      "inv.iter 1 ... chi² = 24.68 (dPhi = 81.24%) lam: 20\n",
      "--------------------------------------------------------------------------------\n",
      "inv.iter 2 ... chi² = 10.36 (dPhi = 54.54%) lam: 20.0\n",
      "--------------------------------------------------------------------------------\n",
      "inv.iter 3 ... chi² = 5.75 (dPhi = 36.32%) lam: 20.0\n",
      "--------------------------------------------------------------------------------\n",
      "inv.iter 4 ... chi² = 4.23 (dPhi = 18.55%) lam: 20.0\n",
      "--------------------------------------------------------------------------------\n",
      "inv.iter 5 ... chi² = 3.82 (dPhi = 5.58%) lam: 20.0\n",
      "--------------------------------------------------------------------------------\n",
      "inv.iter 6 ... chi² = 3.82 (dPhi = -0.02%) lam: 20.0\n",
      "################################################################################\n",
      "#                Abort criterion reached: dPhi = -0.02 (< 2.0%)                #\n",
      "################################################################################\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "20703 [79.18546888337605,...,69.51598653710946]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.estimateError(relativeError=0.05, absoluteUError=50e-6)\n",
    "mgr = ert.Manager(data)\n",
    "mgr.invert(mesh=mesh, zWeight=0.2, verbose=True)\n",
    "# mgr.saveResult()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mesh: Nodes: 4808 Cells: 20703 Boundaries: 5522\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "80dbd439eacd43bf8e074e4c0c12bd35",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Widget(value=\"<iframe src='http://localhost:59357/index.html?ui=P_0x202a2ee5c40_4&reconnect=auto' style='width…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pd = mgr.paraDomain\n",
    "print(pd)\n",
    "pd[\"res\"] = mgr.model\n",
    "pl, _ = pg.show(pd, label=\"res\", style=\"surface\",\n",
    "                filter={\"threshold\": dict(value=20, scalars=\"res\", method=\"upper\")}, \n",
    "                cMap=\"Spectral_r\", cMin=4, cMax=40, logScale=True, hold=True, notebook=True, backend=\"trame\")\n",
    "pv.drawMesh(pl, pd, label=\"res\", style=\"surface\",\n",
    "            filter={\"threshold\": dict(value=6, scalars=\"res\", method=\"lower\")}, \n",
    "            cMap=\"Spectral_r\", cMin=4, cMax=40, logScale=True, hold=True)\n",
    "pl.camera_position = \"xz\"\n",
    "pl.camera.azimuth = +50\n",
    "pl.camera.elevation = 25\n",
    "pl.camera.zoom(1.4)\n",
    "_ = pl.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Inversion result](screenshot.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pg14",
   "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.9.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}