SimpleModel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys, os\n",
    "import subprocess\n",
    "import matplotlib.pyplot as plt\n",
    "from shutil import copyfile\n",
    "# adjust some settings for matplotlib\n",
    "from matplotlib import rcParams\n",
    "import pynoddy\n",
    "# print rcParams\n",
    "# determine path of repository to set paths corretly below\n",
    "#repo_path = os.path.realpath('../..')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'C:\\\\Users\\\\ahinoamp\\\\AppData\\\\Local\\\\Continuum\\\\anaconda3\\\\lib\\\\site-packages\\\\pynoddy\\\\output.py'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "basepynoddyfile = pynoddy.__file__[:-11]+'experiment/__init__.py'\n",
    "# Read in the file\n",
    "with open(basepynoddyfile, 'r') as file :\n",
    "  filedata = file.read()\n",
    "\n",
    "# Replace the target string\n",
    "filedata = filedata.replace('from . import util.sampling as Sample', 'from .util import sampling as Sample')\n",
    "\n",
    "# Write the file out again\n",
    "with open(basepynoddyfile, 'w') as file:\n",
    "  file.write(filedata)\n",
    "  \n",
    "target = pynoddy.__file__[:-11]+'output.py'\n",
    "\n",
    "source = 'output.py'\n",
    "\n",
    "copyfile(source, target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pynoddy.history\n",
    "import pynoddy.experiment\n",
    "import pynoddy.events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'eval' is not recognized as an internal or external command,\n",
      "operable program or batch file.\n"
     ]
    }
   ],
   "source": [
    "folder = !pwd\n",
    "folder = folder[0]\n",
    "noddyEXE = folder+'/noddy.exe'\n",
    "strV = 'chmod 777 '+noddyEXE\n",
    "!eval \"$strV\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cube size is 250.0\n"
     ]
    },
    {
     "ename": "FileNotFoundError",
     "evalue": "[WinError 2] The system cannot find the file specified",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-74217a542cef>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     13\u001b[0m print(subprocess.Popen([noddyEXE, modelfile, output_name, 'ALL'], \n\u001b[0;32m     14\u001b[0m                        \u001b[0mshell\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msubprocess\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mPIPE\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 15\u001b[1;33m                        stdout=subprocess.PIPE).stdout.read())\n\u001b[0m",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\subprocess.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)\u001b[0m\n\u001b[0;32m    773\u001b[0m                                 \u001b[0mc2pread\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mc2pwrite\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    774\u001b[0m                                 \u001b[0merrread\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrwrite\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 775\u001b[1;33m                                 restore_signals, start_new_session)\n\u001b[0m\u001b[0;32m    776\u001b[0m         \u001b[1;32mexcept\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    777\u001b[0m             \u001b[1;31m# Cleanup if the child failed starting.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\subprocess.py\u001b[0m in \u001b[0;36m_execute_child\u001b[1;34m(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, unused_restore_signals, unused_start_new_session)\u001b[0m\n\u001b[0;32m   1176\u001b[0m                                          \u001b[0menv\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1177\u001b[0m                                          \u001b[0mos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfspath\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcwd\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mcwd\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;32mNone\u001b[0m \u001b[1;32melse\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1178\u001b[1;33m                                          startupinfo)\n\u001b[0m\u001b[0;32m   1179\u001b[0m             \u001b[1;32mfinally\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1180\u001b[0m                 \u001b[1;31m# Child is launched. Close the parent's copy of those pipe\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mFileNotFoundError\u001b[0m: [WinError 2] The system cannot find the file specified"
     ]
    }
   ],
   "source": [
    "modelfile = 'ThreeLayerTrenton.his'\n",
    "#pynoddyEXEfile = '/content/DataAndCode/noddy_Mac.exe'\n",
    "output_name = 'noddy_out'\n",
    "\n",
    "H1 = pynoddy.history.NoddyHistory(modelfile)\n",
    "H1.change_cube_size(250)\n",
    "H1.events[3].properties['Dip'] = 40\n",
    "H1.events[4].properties['Dip'] = 40\n",
    "H1.events[5].properties['Dip'] = 40\n",
    "H1.events[5].properties['X'] = 8000\n",
    "H1.write_history(modelfile)\n",
    "print('cube size is ' +str(H1.get_cube_size()))\n",
    "print(subprocess.Popen([noddyEXE, modelfile, output_name, 'ALL'], \n",
    "                       shell=False, stderr=subprocess.PIPE, \n",
    "                       stdout=subprocess.PIPE).stdout.read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " apt.txt\t\t        noddy_Mac\r\n",
      " arr.pvd\t\t        noddy_out.dxf\r\n",
      "'auto_rendering (1).ipynb'      noddy_out.g00\r\n",
      " auto_rendering.ipynb\t        noddy_out.g01\r\n",
      " core.python.140.1578098860     noddy_out.g02\r\n",
      " core.python.168.1578098868     noddy_out.g12\r\n",
      " core.python.53.1578098844      noddy_out.grv\r\n",
      " DisplacementTransferZone.his   noddy_out.mag\r\n",
      " DTZ_3d.png\t\t        noddy_out_play.vtr\r\n",
      " environment.yml\t        noddy_out.vtr\r\n",
      " kkm.vtr\t\t        output.py\r\n",
      " kk.pvtr\t\t        Points.csv\r\n",
      " legacy.vtk\t\t        postBuild\r\n",
      "'Main (1)-Copy1.ipynb'\t        ReallyNiceTrentonModel.his\r\n",
      "'Main (1).ipynb'\t        setup.py\r\n",
      " Main.ipynb\t\t        ThreeLayerTrenton.his\r\n",
      " noddyBatchProgress.txt         WellNamesPaths.csv\r\n",
      " noddy.exe\r\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "463763\n",
      "1391289\n"
     ]
    }
   ],
   "source": [
    "import ezdxf\n",
    "import pandas as pd\n",
    "filename = 'noddy_out.dxf'\n",
    "doc = ezdxf.readfile(filename)\n",
    "\n",
    "def print_entity(e):\n",
    "    print(\"LINE on layer: %s\\n\" % e.dxf.layer)\n",
    "    print(\"start point1: %s\\n\" % e.dxf.vtx0)\n",
    "    print(\"end point2: %s\\n\" % e.dxf.vtx1)\n",
    "    print(\"end point3: %s\\n\" % e.dxf.vtx2)\n",
    "    print(\"end point4: %s\\n\" % e.dxf.vtx3)\n",
    "\n",
    "\n",
    "msp = doc.modelspace()\n",
    "faces = msp.query('3DFACE')\n",
    "num3Dfaces=len(faces)\n",
    "print(num3Dfaces)\n",
    "\n",
    "# %%\n",
    "\n",
    "\n",
    "i=0\n",
    "x = np.zeros((num3Dfaces*3, 1))\n",
    "y = np.zeros((num3Dfaces*3, 1))\n",
    "z = np.zeros((num3Dfaces*3, 1))\n",
    "name = []\n",
    "for e in faces:\n",
    "    x[i] = e.dxf.vtx0[0]\n",
    "    y[i] = e.dxf.vtx0[1]\n",
    "    z[i] = e.dxf.vtx0[2]\n",
    "    name.append(e.dxf.layer)\n",
    "\n",
    "    x[i+1] = e.dxf.vtx1[0]\n",
    "    y[i+1] = e.dxf.vtx1[1]\n",
    "    z[i+1] = e.dxf.vtx1[2]\n",
    "    name.append(e.dxf.layer)\n",
    "\n",
    "    x[i+2] = e.dxf.vtx2[0]\n",
    "    y[i+2] = e.dxf.vtx2[1]\n",
    "    z[i+2] = e.dxf.vtx2[2]\n",
    "    name.append(e.dxf.layer)\n",
    "    i=i+3\n",
    "print(i)\n",
    "\n",
    "# %%\n",
    "#remove row duplicates        \n",
    "AllPoints = pd.DataFrame({'x': x.reshape((-1,)), 'y':y.reshape((-1,)), 'z': z.reshape((-1,)), 'lith': name})\n",
    "AllPoints = AllPoints.drop_duplicates(subset=['x', 'y', 'z'], keep='first')\n",
    "AllPoints.to_csv('Points.csv', index=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ezdxf\n",
    "import pandas as pd\n",
    "filename = 'noddy_out.dxf'\n",
    "doc = ezdxf.readfile(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "463763\n",
      "1391289\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import meshio\n",
    "\n",
    "msp = doc.modelspace()\n",
    "faces = msp.query('3DFACE')\n",
    "num3Dfaces=len(faces)\n",
    "print(num3Dfaces)\n",
    "points = np.zeros((num3Dfaces*3, 3))\n",
    "\n",
    "cells = np.zeros((num3Dfaces, 3),dtype ='int')\n",
    "cell_data=[]\n",
    "\n",
    "i=0\n",
    "faceCounter=0\n",
    "name = []\n",
    "for e in faces:\n",
    "    points[i, 0] = e.dxf.vtx0[0]\n",
    "    points[i, 1] = e.dxf.vtx0[1]\n",
    "    points[i, 2] = e.dxf.vtx0[2]\n",
    "\n",
    "    points[i+1, 0] = e.dxf.vtx1[0]\n",
    "    points[i+1, 1] = e.dxf.vtx1[1]\n",
    "    points[i+1, 2] = e.dxf.vtx1[2]\n",
    "\n",
    "    points[i+2, 0] = e.dxf.vtx2[0]\n",
    "    points[i+2, 1] = e.dxf.vtx2[1]\n",
    "    points[i+2, 2] = e.dxf.vtx2[2]\n",
    "\n",
    "    cells[faceCounter,:]= [i, i+1, i+2]\n",
    "    cell_data.append(e.dxf.layer)\n",
    "    i=i+3\n",
    "    faceCounter=faceCounter+1\n",
    "print(i)\n",
    "\n",
    "cell_data = pd.Series(cell_data)\n",
    "\n",
    "CatCodes = np.zeros((len(cell_data),))\n",
    "filterB = (cell_data.str.contains('B')) \n",
    "filterS = (cell_data.str.contains('S')) \n",
    "\n",
    "CatCodes[filterB]= cell_data.loc[filterB].str[:-16].astype('category').cat.codes\n",
    "CatCodes[filterS]= -1*cell_data.loc[filterS].str[:-12].astype('category').cat.codes\n",
    "\n",
    "#CatCodes[filterB]= 1\n",
    "#CatCodes[filterS]= 2\n",
    "\n",
    "# %%\n",
    "UniqueCodes = np.unique(CatCodes)\n",
    "nCodes = len(UniqueCodes)\n",
    "\n",
    "meshio.write_points_cells(\n",
    "    \"fooColor.vtk\",\n",
    "    points,\n",
    "    cells={'triangle':cells},\n",
    "    cell_data= {'triangle': {'cat':CatCodes}}   \n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "463763\n",
      "1391289\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import meshio\n",
    "\n",
    "msp = doc.modelspace()\n",
    "faces = msp.query('3DFACE')\n",
    "num3Dfaces=len(faces)\n",
    "print(num3Dfaces)\n",
    "points = np.zeros((num3Dfaces*3, 3))\n",
    "\n",
    "cells = np.zeros((num3Dfaces, 3),dtype ='int')\n",
    "cell_data=[]\n",
    "\n",
    "i=0\n",
    "faceCounter=0\n",
    "name = []\n",
    "for e in faces:\n",
    "    points[i, 0] = e.dxf.vtx0[0]\n",
    "    points[i, 1] = e.dxf.vtx0[1]\n",
    "    points[i, 2] = e.dxf.vtx0[2]\n",
    "\n",
    "    points[i+1, 0] = e.dxf.vtx1[0]\n",
    "    points[i+1, 1] = e.dxf.vtx1[1]\n",
    "    points[i+1, 2] = e.dxf.vtx1[2]\n",
    "\n",
    "    points[i+2, 0] = e.dxf.vtx2[0]\n",
    "    points[i+2, 1] = e.dxf.vtx2[1]\n",
    "    points[i+2, 2] = e.dxf.vtx2[2]\n",
    "\n",
    "    cells[faceCounter,:]= [i, i+1, i+2]\n",
    "    cell_data.append(e.dxf.layer)\n",
    "    i=i+3\n",
    "    faceCounter=faceCounter+1\n",
    "print(i)\n",
    "\n",
    "cell_data = pd.Series(cell_data)\n",
    "\n",
    "CatCodes = np.zeros((len(cell_data),))\n",
    "filterB = (cell_data.str.contains('B')) \n",
    "filterS = (cell_data.str.contains('S')) \n",
    "\n",
    "CatCodes[filterB]= cell_data.loc[filterB].str[:-16].astype('category').cat.codes\n",
    "CatCodes[filterS]= -1*cell_data.loc[filterS].str[:-12].astype('category').cat.codes\n",
    "\n",
    "#CatCodes[filterB]= 1\n",
    "#CatCodes[filterS]= 2\n",
    "\n",
    "# %%\n",
    "UniqueCodes = np.unique(CatCodes)\n",
    "nCodes = len(UniqueCodes)\n",
    "\n",
    "UniqueCodes = np.unique(CatCodes)\n",
    "nCodes = len(UniqueCodes)\n",
    "for i in range(nCodes):\n",
    "    filterPoints = CatCodes==UniqueCodes[i]\n",
    "    nCells = np.sum(filterPoints)\n",
    "    Cells_i = np.zeros((nCells, 3),dtype ='int')\n",
    "    cntr = 0\n",
    "    for j in range(nCells):\n",
    "        Cells_i[j]=[cntr, cntr+1, cntr+2]\n",
    "        cntr=cntr+3\n",
    "        \n",
    "    meshio.write_points_cells(\n",
    "        \"fooTry\"+str(i)+\".vtk\",\n",
    "        points[np.repeat(filterPoints,3), :],\n",
    "        cells={'triangle':Cells_i}\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9, 10, 11])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0\n",
    "faceCounter=0\n",
    "name = []\n",
    "for e in faces:\n",
    "    points[i, 0] = e.dxf.vtx0[0]\n",
    "    points[i, 1] = e.dxf.vtx0[1]\n",
    "    points[i, 2] = e.dxf.vtx0[2]\n",
    "\n",
    "    points[i+1, 0] = e.dxf.vtx1[0]\n",
    "    points[i+1, 1] = e.dxf.vtx1[1]\n",
    "    points[i+1, 2] = e.dxf.vtx1[2]\n",
    "\n",
    "    points[i+2, 0] = e.dxf.vtx2[0]\n",
    "    points[i+2, 1] = e.dxf.vtx2[1]\n",
    "    points[i+2, 2] = e.dxf.vtx2[2]\n",
    "\n",
    "    cells[faceCounter,:]= [i, i+1, i+2]\n",
    "    cell_data.append(e.dxf.layer)\n",
    "    i=i+3\n",
    "    faceCounter=faceCounter+1\n",
    "print(i)\n",
    "\n",
    "cell_data = pd.Series(cell_data)\n",
    "\n",
    "CatCodes = np.zeros((len(cell_data),))\n",
    "filterB = (cell_data.str.contains('B')) \n",
    "filterS = (cell_data.str.contains('S')) \n",
    "\n",
    "CatCodes[filterB]= cell_data.loc[filterB].str[:-16].astype('category').cat.codes\n",
    "CatCodes[filterS]= -1*cell_data.loc[filterS].str[:-12].astype('category').cat.codes\n",
    "\n",
    "CatCodes[filterB]= 1\n",
    "CatCodes[filterS]= 2\n",
    "\n",
    "# %%\n",
    "UniqueCodes = np.unique(CatCodes)\n",
    "nCodes = len(UniqueCodes)\n",
    "for i in range(nCodes):\n",
    "    filterPoints = CatCodes==UniqueCodes[i]\n",
    "    meshio.write_points_cells(\n",
    "        \"fooTry\"+str(i)+\".vtk\",\n",
    "        points[np.repeat(filterPoints,3), :],\n",
    "        cells={'triangle':cells[filterPoints, :]}\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "463763"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "UniqueCodes = np.unique(CatCodes)\n",
    "nCodes = len(UniqueCodes)\n",
    "for i in range(nCodes):\n",
    "    filterPoints = CatCodes==UniqueCodes[i]\n",
    "    meshio.write_points_cells(\n",
    "        \"foo\"+str(i)+\".vtk\",\n",
    "        points[np.repeat(filterPoints,3), :],\n",
    "        cells={'triangle':cells[filterPoints, :]}\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1391289"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(np.repeat(filterPoints,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-49-fd9c056d4157>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mfilterS\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mcell_data\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcontains\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'S'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m \u001b[0mCatCodes\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mfilterB\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m \u001b[0mcell_data\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mfilterB\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m16\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'category'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcat\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcodes\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      8\u001b[0m \u001b[0mCatCodes\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mfilterS\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mcell_data\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mfilterS\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m12\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'category'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcat\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcodes\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mIndexError\u001b[0m: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices"
     ]
    }
   ],
   "source": [
    "cell_data = pd.Series(cell_data)\n",
    "\n",
    "CatCodes = np.zeros((len(cell_data),))\n",
    "filterB = (cell_data.str.contains('B')) \n",
    "filterS = (cell_data.str.contains('S')) \n",
    "\n",
    "CatCodes[filterB]= cell_data[filterB].str[:-16].astype('category').cat.codes\n",
    "CatCodes[filterS]= -1*cell_data[filterS].str[:-12].astype('category').cat.codes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lith    lith    lith    lith    0         S_00_02_001_...\n",
       "dtype: object"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cell_data[0:1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_name = 'noddy_out'\n",
    "N1 = pynoddy.output.NoddyOutput(output_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The model has an extent of 14000 m in x-direction, with 56 cells of width 250 m\n"
     ]
    }
   ],
   "source": [
    "print(\"The model has an extent of %.0f m in x-direction, with %d cells of width %.0f m\" %\n",
    "      (N1.extent_x, N1.nx, N1.delx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "N1.export_to_vtk()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polyfit2d(x, y, z, order=3):\n",
    "    ncols = (order + 1)**2\n",
    "    G = np.zeros((x.size, ncols))\n",
    "    ij = itertools.product(range(order+1), range(order+1))\n",
    "#    ij = [tup for tup in ij if sum(tup) <= order]\n",
    "    for k, (i,j) in enumerate(ij):\n",
    "        G[:,k] = x**i * y**j\n",
    "    m, _, _, _ = np.linalg.lstsq(G, z)\n",
    "    return m\n",
    "\n",
    "def polyval2d(x, y, m):\n",
    "    order = int(np.sqrt(len(m))) - 1\n",
    "    ij = itertools.product(range(order+1), range(order+1))\n",
    "    z = np.zeros_like(x)\n",
    "    for a, (i,j) in zip(m, ij):\n",
    "        z += a * x**i * y**j\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(23275, 3)\n",
      "Recalculating normals for 4931 Points, sample size= 98\n",
      "Calculating mesh from points with R = 542.9252451182724\n",
      "(23275, 3)\n",
      "Recalculating normals for 4517 Points, sample size= 90\n",
      "Calculating mesh from points with R = 555.1824375865258\n",
      "(23275, 3)\n",
      "Recalculating normals for 2422 Points, sample size= 48\n",
      "Calculating mesh from points with R = 348.9759449929405\n",
      "(23275, 3)\n",
      "Recalculating normals for 2510 Points, sample size= 50\n",
      "Calculating mesh from points with R = 348.9759449929405\n",
      "(23275, 3)\n",
      "Recalculating normals for 2607 Points, sample size= 52\n",
      "Calculating mesh from points with R = 351.4607599877364\n",
      "(23275, 3)\n",
      "Recalculating normals for 2370 Points, sample size= 47\n",
      "Calculating mesh from points with R = 348.9759449929405\n",
      "(23275, 3)\n",
      "Recalculating normals for 1505 Points, sample size= 30\n",
      "Calculating mesh from points with R = 301.82733317800216\n",
      "(23275, 3)\n",
      "Recalculating normals for 2413 Points, sample size= 48\n",
      "Calculating mesh from points with R = 343.9178056904336\n",
      "\u001b[1m\u001b[31m❌ pointColors Error: nr. of scalars != nr. of points 7712 7714\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from vtkplotter import *\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import itertools\n",
    "import k3d\n",
    "from scipy.interpolate import griddata\n",
    "import vtk\n",
    "\n",
    "def GetCat(strCat):\n",
    "    CatCodes = np.zeros((len(strCat),))\n",
    "    filterB = (strCat['lith'].str.contains('B')) \n",
    "    filterS = (strCat['lith'].str.contains('S')) \n",
    "    \n",
    "    CatCodes[filterB]= strCat.loc[filterB, 'lith'].str[:-16].astype('category').cat.codes\n",
    "    CatCodes[filterS]= -1*strCat.loc[filterS, 'lith'].str[:-12].astype('category').cat.codes\n",
    "\n",
    "#    CatCodes[filterB]= 0\n",
    "#    CatCodes[filterS]= 1\n",
    "\n",
    "    return CatCodes\n",
    "\n",
    "\n",
    "embedWindow('k3d') #uncomment to pop the usual vtk rendering window, press Esc to close it.\n",
    "\n",
    "positions = np.random.randn(5000,3)\n",
    "\n",
    "data = pd.read_csv('Points.csv')\n",
    "data['Cat'] = GetCat(data)\n",
    "\n",
    "positions2 = data[['x', 'y','z']].values\n",
    "scalars = data['Cat']\n",
    "\n",
    "plot = Plotter()\n",
    "#plot = k3d.plot()\n",
    "\n",
    "plot.xtitle = '\\sum_i x_i'\n",
    "\n",
    "\n",
    "listSurfaces = np.unique(data['Cat'])\n",
    "nSurfaces = len(listSurfaces)\n",
    "colors = pl.cm.jet(np.linspace(0,1,nSurfaces))\n",
    "\n",
    "for i in range(nSurfaces):\n",
    "    surfaceId = listSurfaces[i]\n",
    "    filterSurface = data['Cat']==surfaceId\n",
    "    print(np.shape(positions2))\n",
    "\n",
    "    points = Points(positions2[filterSurface]).color(colors[i][0:3])\n",
    "    # color vertices based on their scalar value with any matplotlib color map\n",
    "#    points.pointColor(colors[i])\n",
    "#    plot += points\n",
    "\n",
    "    reco = recoSurface(points, bins=128).legend(\"surf reco\")\n",
    "#    plt_surface = k3d.surface(zz, xmin=x.min(), xmax=x.max(), ymin=y.min(), ymax=y.max())\n",
    "\n",
    "    plot += reco\n",
    "\n",
    "filename = 'kkm.vtr'\n",
    "e=load(filename)\n",
    "#e = loadRectilinearGrid(filename)\n",
    "#reader = vtk.vtkRectilinearGridReader()\n",
    "#reader.SetFileName(filename)\n",
    "#reader.Update()\n",
    "#routput = reader.GetOutput()\n",
    "#e = Actor(routput)\n",
    "\n",
    "value = e.info.keys()\n",
    "scals = e.scalars('geology')\n",
    "e.pointColors(scals, cmap='Set3')\n",
    "\n",
    "plot += e\n",
    "\n",
    "#plot.show()\n",
    "#plot.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unique(data['Cat'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import k3d\n",
    "\n",
    "plot = k3d.plot()\n",
    "\n",
    "Nx, Ny = 40, 80\n",
    "xmin, xmax, ymin, ymax = -3, 3, 0, 3\n",
    "\n",
    "x = np.linspace(xmin, xmax, Nx, dtype=np.float32)\n",
    "y = np.linspace(ymin, ymax, Ny, dtype=np.float32)\n",
    "x, y = np.meshgrid(x, y)\n",
    "\n",
    "heights = np.sin(x**2 + y**2)\n",
    "\n",
    "plt_surface = k3d.surface(heights, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)\n",
    "plot += plt_surface\n",
    "\n",
    "plot.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[31m❌ pointColors Error: nr. of scalars != nr. of points 7712 7714\u001b[0m\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'plot' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-ee90293ec0aa>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpointColors\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Set3'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mplot\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mplot\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'plot' is not defined"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from vtkplotter import *\n",
    "\n",
    "filename = 'noddy_out.vtr'\n",
    "e=load(filename)\n",
    "scals = e.scalars('geology')\n",
    "e.pointColors(scals, cmap='Set3')\n",
    "\n",
    "plot += e\n",
    "\n",
    "plot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(40320, 3)\n",
      "(20, 36, 56)\n",
      "(40320, 3)\n",
      "(40320,)\n",
      "(40, 40)\n",
      "(40, 40)\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "SetDimensions argument 1: expected a sequence of 3 values, got 2 values",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-119-6eec1b7355ef>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     45\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     46\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 47\u001b[0;31m \u001b[0mvol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mVolume\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnewc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'prism'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     48\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[0mlego\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvol\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlegosurface\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/srv/conda/envs/notebook/lib/python3.7/site-packages/vtkplotter/actors.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, inputobj, c, alpha, alphaGradient, mode, origin, spacing, mapperType)\u001b[0m\n\u001b[1;32m   3444\u001b[0m             \u001b[0mvarr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSetName\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'input_scalars'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3445\u001b[0m             \u001b[0mimg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvtk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvtkImageData\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3446\u001b[0;31m             \u001b[0mimg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSetDimensions\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3447\u001b[0m             \u001b[0mimg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGetPointData\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSetScalars\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvarr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3448\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0;34m\"ImageData\"\u001b[0m \u001b[0;32min\u001b[0m \u001b[0minputtype\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: SetDimensions argument 1: expected a sequence of 3 values, got 2 values"
     ]
    }
   ],
   "source": [
    "from vtkplotter import *\n",
    "\n",
    "embedWindow('k3d') #uncomment to pop the usual vtk rendering window, press Esc to close it.\n",
    "\n",
    "N1 = pynoddy.output.NoddyOutput(output_name)\n",
    "\n",
    "x = np.arange(0, N1.extent_x, N1.delx, dtype='float64')\n",
    "y = np.arange(0, N1.extent_y, N1.dely, dtype='float64')\n",
    "z = np.arange(0, N1.extent_z, N1.delz, dtype='float64')\n",
    "color = N1.block\n",
    "\n",
    "#PointsPlot = np.concatenate((x.reshape(-1,1), y.reshape(-1,1), z.reshape(-1,1)), axis=1)\n",
    "AAA= np.meshgrid(x,y,z)\n",
    "\n",
    "x=AAA[0].reshape(-1,1)\n",
    "y=AAA[1].reshape(-1,1)\n",
    "z=AAA[2].reshape(-1,1)\n",
    "\n",
    "Coord = np.concatenate((x,y,z), axis=1)\n",
    "\n",
    "plot = Plotter()\n",
    "print(np.shape(Coord))\n",
    "points = Points(Coord).pointColors(np.transpose(color, axes =(1, 0, 2)).reshape(-1,1), cmap='Set3')\n",
    "#plot += points\n",
    "\n",
    "print(np.shape(np.transpose(color, axes =(2, 1, 0))))\n",
    "             \n",
    "#vol = Volume(np.transpose(color, axes =(2, 1, 0)), c='prism', spacing=[N1.delx, N1.dely,N1.delz])\n",
    "\n",
    "xf = np.linspace(x.min(), x.max(), 40, dtype=np.float32)\n",
    "yf = np.linspace(y.min(), y.max(), 40, dtype=np.float32)\n",
    "xx, yy = np.meshgrid(xf, yf)\n",
    "\n",
    "##    zfun_smooth_rbf = interp.Rbf(x, y, z, function='cubic', smooth=0)  # default smooth=0 for interpolation\n",
    "##    zz = zfun_smooth_rbf(xx, yy)  # not really a function, but a callable class instance\n",
    "\n",
    "pointsE=np.concatenate((x.reshape(-1,1), y.reshape(-1,1), z.reshape(-1,1)), axis=1)\n",
    "print(np.shape(pointsE))\n",
    "print(np.shape(color.reshape(-1,)))\n",
    "print(np.shape(xx))\n",
    "print(np.shape(yy))\n",
    "\n",
    "newc = griddata(pointsE, color.reshape(-1,), (xx, yy, zz), method='nearest')\n",
    "\n",
    "\n",
    "\n",
    "vol = Volume(newc, c='prism')\n",
    "\n",
    "lego = vol.legosurface(-1, 4)\n",
    "\n",
    "\n",
    "plot += lego\n",
    "\n",
    "print('scalar min, max =', np.min(color), np.max(color))\n",
    "\n",
    "plot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[[    0.,     0.,     0., ...,     0.,     0.,     0.],\n",
       "         [  250.,   250.,   250., ...,   250.,   250.,   250.],\n",
       "         [  500.,   500.,   500., ...,   500.,   500.,   500.],\n",
       "         ...,\n",
       "         [13500., 13500., 13500., ..., 13500., 13500., 13500.],\n",
       "         [13750., 13750., 13750., ..., 13750., 13750., 13750.],\n",
       "         [14000., 14000., 14000., ..., 14000., 14000., 14000.]],\n",
       " \n",
       "        [[    0.,     0.,     0., ...,     0.,     0.,     0.],\n",
       "         [  250.,   250.,   250., ...,   250.,   250.,   250.],\n",
       "         [  500.,   500.,   500., ...,   500.,   500.,   500.],\n",
       "         ...,\n",
       "         [13500., 13500., 13500., ..., 13500., 13500., 13500.],\n",
       "         [13750., 13750., 13750., ..., 13750., 13750., 13750.],\n",
       "         [14000., 14000., 14000., ..., 14000., 14000., 14000.]],\n",
       " \n",
       "        [[    0.,     0.,     0., ...,     0.,     0.,     0.],\n",
       "         [  250.,   250.,   250., ...,   250.,   250.,   250.],\n",
       "         [  500.,   500.,   500., ...,   500.,   500.,   500.],\n",
       "         ...,\n",
       "         [13500., 13500., 13500., ..., 13500., 13500., 13500.],\n",
       "         [13750., 13750., 13750., ..., 13750., 13750., 13750.],\n",
       "         [14000., 14000., 14000., ..., 14000., 14000., 14000.]],\n",
       " \n",
       "        ...,\n",
       " \n",
       "        [[    0.,     0.,     0., ...,     0.,     0.,     0.],\n",
       "         [  250.,   250.,   250., ...,   250.,   250.,   250.],\n",
       "         [  500.,   500.,   500., ...,   500.,   500.,   500.],\n",
       "         ...,\n",
       "         [13500., 13500., 13500., ..., 13500., 13500., 13500.],\n",
       "         [13750., 13750., 13750., ..., 13750., 13750., 13750.],\n",
       "         [14000., 14000., 14000., ..., 14000., 14000., 14000.]],\n",
       " \n",
       "        [[    0.,     0.,     0., ...,     0.,     0.,     0.],\n",
       "         [  250.,   250.,   250., ...,   250.,   250.,   250.],\n",
       "         [  500.,   500.,   500., ...,   500.,   500.,   500.],\n",
       "         ...,\n",
       "         [13500., 13500., 13500., ..., 13500., 13500., 13500.],\n",
       "         [13750., 13750., 13750., ..., 13750., 13750., 13750.],\n",
       "         [14000., 14000., 14000., ..., 14000., 14000., 14000.]],\n",
       " \n",
       "        [[    0.,     0.,     0., ...,     0.,     0.,     0.],\n",
       "         [  250.,   250.,   250., ...,   250.,   250.,   250.],\n",
       "         [  500.,   500.,   500., ...,   500.,   500.,   500.],\n",
       "         ...,\n",
       "         [13500., 13500., 13500., ..., 13500., 13500., 13500.],\n",
       "         [13750., 13750., 13750., ..., 13750., 13750., 13750.],\n",
       "         [14000., 14000., 14000., ..., 14000., 14000., 14000.]]]),\n",
       " array([[[   0.,    0.,    0., ...,    0.,    0.,    0.],\n",
       "         [   0.,    0.,    0., ...,    0.,    0.,    0.],\n",
       "         [   0.,    0.,    0., ...,    0.,    0.,    0.],\n",
       "         ...,\n",
       "         [   0.,    0.,    0., ...,    0.,    0.,    0.],\n",
       "         [   0.,    0.,    0., ...,    0.,    0.,    0.],\n",
       "         [   0.,    0.,    0., ...,    0.,    0.,    0.]],\n",
       " \n",
       "        [[ 250.,  250.,  250., ...,  250.,  250.,  250.],\n",
       "         [ 250.,  250.,  250., ...,  250.,  250.,  250.],\n",
       "         [ 250.,  250.,  250., ...,  250.,  250.,  250.],\n",
       "         ...,\n",
       "         [ 250.,  250.,  250., ...,  250.,  250.,  250.],\n",
       "         [ 250.,  250.,  250., ...,  250.,  250.,  250.],\n",
       "         [ 250.,  250.,  250., ...,  250.,  250.,  250.]],\n",
       " \n",
       "        [[ 500.,  500.,  500., ...,  500.,  500.,  500.],\n",
       "         [ 500.,  500.,  500., ...,  500.,  500.,  500.],\n",
       "         [ 500.,  500.,  500., ...,  500.,  500.,  500.],\n",
       "         ...,\n",
       "         [ 500.,  500.,  500., ...,  500.,  500.,  500.],\n",
       "         [ 500.,  500.,  500., ...,  500.,  500.,  500.],\n",
       "         [ 500.,  500.,  500., ...,  500.,  500.,  500.]],\n",
       " \n",
       "        ...,\n",
       " \n",
       "        [[8500., 8500., 8500., ..., 8500., 8500., 8500.],\n",
       "         [8500., 8500., 8500., ..., 8500., 8500., 8500.],\n",
       "         [8500., 8500., 8500., ..., 8500., 8500., 8500.],\n",
       "         ...,\n",
       "         [8500., 8500., 8500., ..., 8500., 8500., 8500.],\n",
       "         [8500., 8500., 8500., ..., 8500., 8500., 8500.],\n",
       "         [8500., 8500., 8500., ..., 8500., 8500., 8500.]],\n",
       " \n",
       "        [[8750., 8750., 8750., ..., 8750., 8750., 8750.],\n",
       "         [8750., 8750., 8750., ..., 8750., 8750., 8750.],\n",
       "         [8750., 8750., 8750., ..., 8750., 8750., 8750.],\n",
       "         ...,\n",
       "         [8750., 8750., 8750., ..., 8750., 8750., 8750.],\n",
       "         [8750., 8750., 8750., ..., 8750., 8750., 8750.],\n",
       "         [8750., 8750., 8750., ..., 8750., 8750., 8750.]],\n",
       " \n",
       "        [[9000., 9000., 9000., ..., 9000., 9000., 9000.],\n",
       "         [9000., 9000., 9000., ..., 9000., 9000., 9000.],\n",
       "         [9000., 9000., 9000., ..., 9000., 9000., 9000.],\n",
       "         ...,\n",
       "         [9000., 9000., 9000., ..., 9000., 9000., 9000.],\n",
       "         [9000., 9000., 9000., ..., 9000., 9000., 9000.],\n",
       "         [9000., 9000., 9000., ..., 9000., 9000., 9000.]]]),\n",
       " array([[[   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         ...,\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.]],\n",
       " \n",
       "        [[   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         ...,\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.]],\n",
       " \n",
       "        [[   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         ...,\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.]],\n",
       " \n",
       "        ...,\n",
       " \n",
       "        [[   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         ...,\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.]],\n",
       " \n",
       "        [[   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         ...,\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.]],\n",
       " \n",
       "        [[   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         ...,\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.],\n",
       "         [   0.,  250.,  500., ..., 4500., 4750., 5000.]]])]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PointsPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(23275, 3)\n",
      "(4931, 2)\n",
      "(4931,)\n",
      "(40, 40)\n",
      "(40, 40)\n",
      "(23275, 3)\n",
      "(4517, 2)\n",
      "(4517,)\n",
      "(40, 40)\n",
      "(40, 40)\n",
      "(23275, 3)\n",
      "(2422, 2)\n",
      "(2422,)\n",
      "(40, 40)\n",
      "(40, 40)\n",
      "(23275, 3)\n",
      "(2510, 2)\n",
      "(2510,)\n",
      "(40, 40)\n",
      "(40, 40)\n",
      "(23275, 3)\n",
      "(2607, 2)\n",
      "(2607,)\n",
      "(40, 40)\n",
      "(40, 40)\n",
      "(23275, 3)\n",
      "(2370, 2)\n",
      "(2370,)\n",
      "(40, 40)\n",
      "(40, 40)\n",
      "(23275, 3)\n",
      "(1505, 2)\n",
      "(1505,)\n",
      "(40, 40)\n",
      "(40, 40)\n",
      "(23275, 3)\n",
      "(2413, 2)\n",
      "(2413,)\n",
      "(40, 40)\n",
      "(40, 40)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0318e9873686418fbbc157f4c7f314c2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#def voxels(voxels, color_map=nice_colors, wireframe=False, outlines=True, outlines_color=0, opacity=1.0, name=None,\n",
    "#           compression_level=0, bounds=None,\n",
    "#           **kwargs):\n",
    "\n",
    "import numpy as np\n",
    "import k3d\n",
    "import scipy.interpolate as interp\n",
    "\n",
    "embedWindow('itk') #uncomment to pop the usual vtk rendering window, press Esc to close it.\n",
    "\n",
    "\n",
    "plot = k3d.plot()\n",
    "\n",
    "\n",
    "voxels = k3d.voxels(np.transpose(color, axes =(2, 1, 0)), bounds=[0, 12000, 0, 7000, 0, 5000])\n",
    "plot += voxels\n",
    "\n",
    "listSurfaces = np.unique(data['Cat'])\n",
    "nSurfaces = len(listSurfaces)\n",
    "colors = pl.cm.jet(np.linspace(0,1,nSurfaces))\n",
    "\n",
    "for i in range(nSurfaces):\n",
    "    surfaceId = listSurfaces[i]\n",
    "    filterSurface = data['Cat']==surfaceId\n",
    "    print(np.shape(positions2))\n",
    "\n",
    "    points = Points(positions2[filterSurface]).color(colors[i][0:3])\n",
    "    # color vertices based on their scalar value with any matplotlib color map\n",
    "#    points.pointColor(colors[i])\n",
    "#    plot += points\n",
    "    x = positions2[filterSurface, 0]\n",
    "    y = positions2[filterSurface, 1]\n",
    "    z = positions2[filterSurface, 2]\n",
    "\n",
    "    xf = np.linspace(x.min(), x.max(), 40, dtype=np.float32)\n",
    "    yf = np.linspace(y.min(), y.max(), 40, dtype=np.float32)\n",
    "    xx, yy = np.meshgrid(xf, yf)\n",
    "    \n",
    "##    zfun_smooth_rbf = interp.Rbf(x, y, z, function='cubic', smooth=0)  # default smooth=0 for interpolation\n",
    "##    zz = zfun_smooth_rbf(xx, yy)  # not really a function, but a callable class instance\n",
    "\n",
    "    zz = griddata(np.concatenate((x.reshape(-1,1), y.reshape(-1,1)), axis=1), z, (xx, yy), method='nearest')\n",
    "\n",
    "    print(np.shape(np.concatenate((x.reshape(-1,1), y.reshape(-1,1)), axis=1)))\n",
    "    print(np.shape(z))\n",
    "    print(np.shape(xx))\n",
    "    print(np.shape(yy))\n",
    "\n",
    "    #    reco = recoSurface(points, bins=128).legend(\"surf reco\")\n",
    "    plt_surface = k3d.surface(zz, xmin=x.min(), xmax=x.max(), ymin=y.min(), ymax=y.max())\n",
    "\n",
    "    plot += plt_surface\n",
    "\n",
    "\n",
    "plot.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14000.0"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N1.delx*56"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scalar min, max = 0.0 1.0\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "00dfeecf7bf14d8a8d7a6b0092a1356f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(cmap='jet', geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_point…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from __future__ import print_function, division\n",
    "import numpy as np\n",
    "from vtkplotter import *\n",
    "embedWindow('itkwidgets') # backends are: itkwidgets, k3d or False\n",
    "\n",
    "X, Y, Z = np.mgrid[:30, :30, :30]\n",
    "# scaled distance from the center at (15, 15, 15)\n",
    "scalar_field = ((X-15)**2 + (Y-15)**2 + (Z-15)**2)/225/3\n",
    "print('scalar min, max =', np.min(scalar_field), np.max(scalar_field))\n",
    "\n",
    "vol = Volume(scalar_field)\n",
    "lego = vol.legosurface(vmin=0, vmax=1.1)\n",
    "\n",
    "#show(vol.c('magma'))\n",
    "show(lego)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import meshio\n",
    "import numpy as np\n",
    "\n",
    "points = np.array([\n",
    "    [0.0, 0.0, 0.0],\n",
    "    [0.0, 1.0, 0.0],\n",
    "    [0.0, 0.0, 1.0],\n",
    "    ])\n",
    "cells = {\n",
    "    \"triangle\": np.array([\n",
    "        [0, 1, 2]\n",
    "        ])\n",
    "    }\n",
    "meshio.write_points_cells(\n",
    "    \"foo.vtk\",\n",
    "    points,\n",
    "    cells,\n",
    "    # Optionally provide extra data on points, cells, etc.\n",
    "    # point_data=point_data,\n",
    "    # cell_data=cell_data,\n",
    "    # field_data=field_data\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "points = np.array([\n",
    "    [0.0, 0.0, 0.0],\n",
    "[1.0, 0.0, 0.0],\n",
    "[2.0, 0.0, 0.0],\n",
    "[0.0, 1.0, 0.0],\n",
    "[1.0, 1.0, 0.0],\n",
    "[2.0, 1.0, 2]])\n",
    "\n",
    "# Define connectivity or vertices that belongs to each element\n",
    "cells = {\n",
    "    \"triangle\": np.array([\n",
    "        [0, 1, 3]\n",
    "        , [1, 4, 3]\n",
    "        ]),\n",
    "    \"quad\": np.array([\n",
    "        [1, 2, 5, 4]\n",
    "        ])\n",
    "    }\n",
    "\n",
    "meshio.write_points_cells(\n",
    "    \"foo.vtk\",\n",
    "    points,\n",
    "    cells,\n",
    "    # Optionally provide extra data on points, cells, etc.\n",
    "    # point_data=point_data,\n",
    "    # cell_data=cell_data,\n",
    "    # field_data=field_data\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'triangle': array([[0, 1, 3],\n",
       "        [1, 4, 3]]), 'quad': array([[1, 2, 5, 4]])}"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from vtkplotter import *\n",
    "import numpy as np\n",
    "#import matplotlib.pylab as pl\n",
    "import itertools\n",
    "import k3d\n",
    "from scipy.interpolate import griddata\n",
    "import vtk\n",
    "\n",
    "plot = Plotter()\n",
    "\n",
    "embedWindow('k3d')\n",
    "\n",
    "filename = 'foo28.vtk'\n",
    "e=load(filename)\n",
    "\n",
    "plot += e\n",
    "\n",
    "plot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d72e77b1380f454ebf76a3a7501f6e0f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Plot(antialias=3, axes=['x', 'y', 'z'], axes_helper=1.0, background_color=16777215, camera=[6146.346336284672,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from vtkplotter import *\n",
    "import numpy as np\n",
    "#import matplotlib.pylab as pl\n",
    "import itertools\n",
    "import k3d\n",
    "from scipy.interpolate import griddata\n",
    "import vtk\n",
    "import matplotlib\n",
    "import matplotlib.pylab as pl\n",
    "\n",
    "plot = Plotter()\n",
    "\n",
    "embedWindow('k3d')\n",
    "\n",
    "nSurfaces = 29\n",
    "colors = pl.cm.jet(np.linspace(0,1,nSurfaces))\n",
    "\n",
    "for i in range(nSurfaces):\n",
    "    print(i)\n",
    "    filename = 'fooTry'+str(i)+'.vtk'\n",
    "    e=load(filename).c(colors[i, 0:3])\n",
    "\n",
    "    plot += e\n",
    "\n",
    "plot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: matplotlib in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (3.1.1)\n",
      "Requirement already satisfied: cycler>=0.10 in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (from matplotlib) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (from matplotlib) (1.1.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (from matplotlib) (2.4.2)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (from matplotlib) (2.8.0)\n",
      "Requirement already satisfied: numpy>=1.11 in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (from matplotlib) (1.16.5)\n",
      "Requirement already satisfied: six in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (from cycler>=0.10->matplotlib) (1.12.0)\n",
      "Requirement already satisfied: setuptools in c:\\users\\ahinoamp\\appdata\\local\\continuum\\anaconda3\\lib\\site-packages (from kiwisolver>=1.0.1->matplotlib) (41.4.0)\n"
     ]
    }
   ],
   "source": [
    "!pip install matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'matplotlib'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-32-0484cd13f94d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'matplotlib'"
     ]
    }
   ],
   "source": [
    "import matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "70d5b35d8ca944c3a06524cc7da2eab0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Plot(antialias=3, axes=['x', 'y', 'z'], axes_helper=1.0, background_color=16777215, camera=[6146.346336284672,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from vtkplotter import *\n",
    "import numpy as np\n",
    "#import matplotlib.pylab as pl\n",
    "import itertools\n",
    "import k3d\n",
    "from scipy.interpolate import griddata\n",
    "import vtk\n",
    "\n",
    "plot = Plotter()\n",
    "\n",
    "embedWindow('k3d')\n",
    "i=1\n",
    "print(i)\n",
    "filename = 'fooColor.vtk'\n",
    "e=load(filename)\n",
    "\n",
    "plot += e\n",
    "\n",
    "plot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}