Main.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "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": [
    {
     "ename": "NameError",
     "evalue": "name 'pynoddy' 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-3-ef23292a45e1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mbasepynoddyfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpynoddy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__file__\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m11\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'experiment/__init__.py'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;31m# Read in the file\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbasepynoddyfile\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mfile\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m   \u001b[0mfiledata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfile\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\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[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'pynoddy' is not defined"
     ]
    }
   ],
   "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": 5,
   "metadata": {},
   "outputs": [],
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cube size is 250.0\n",
      "b''\n"
     ]
    }
   ],
   "source": [
    "modelfile = 'ReallyNiceTrentonModel.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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "apt.txt\t\t\t      noddy_Mac      noddy_out.vtr\r\n",
      "auto_rendering.ipynb\t      noddy_out.dxf  output.py\r\n",
      "core.python.23.1578085425     noddy_out.g00  Points.csv\r\n",
      "DisplacementTransferZone.his  noddy_out.g01  postBuild\r\n",
      "DTZ_3d.png\t\t      noddy_out.g02  ReallyNiceTrentonModel.his\r\n",
      "environment.yml\t\t      noddy_out.g12  setup.py\r\n",
      "noddyBatchProgress.txt\t      noddy_out.grv  WellNamesPaths.csv\r\n",
      "noddy.exe\t\t      noddy_out.mag\r\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "170025\n",
      "510075\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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_name = 'noddy_out'\n",
    "N1 = pynoddy.output.NoddyOutput(output_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "N1.export_to_vtk()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Recalculating normals for 4931 Points, sample size= 98\n",
      "Calculating mesh from points with R = 542.9252451182724\n",
      "Recalculating normals for 4858 Points, sample size= 97\n",
      "Calculating mesh from points with R = 542.9057500487331\n",
      "Recalculating normals for 4814 Points, sample size= 96\n",
      "Calculating mesh from points with R = 542.8936039204636\n",
      "Recalculating normals for 4826 Points, sample size= 96\n",
      "Calculating mesh from points with R = 542.8735178764671\n",
      "Recalculating normals for 4025 Points, sample size= 80\n",
      "Calculating mesh from points with R = 542.7531031516357\n",
      "Recalculating normals for 4552 Points, sample size= 91\n",
      "Calculating mesh from points with R = 542.7467424395546\n",
      "Recalculating normals for 4615 Points, sample size= 92\n",
      "Calculating mesh from points with R = 542.8352232599074\n",
      "Recalculating normals for 4435 Points, sample size= 88\n",
      "Calculating mesh from points with R = 542.8148496533743\n",
      "Recalculating normals for 4446 Points, sample size= 88\n",
      "Calculating mesh from points with R = 542.8053952966469\n",
      "Recalculating normals for 2769 Points, sample size= 55\n",
      "Calculating mesh from points with R = 542.7121467033894\n",
      "Recalculating normals for 4405 Points, sample size= 88\n",
      "Calculating mesh from points with R = 542.7074191413399\n",
      "Recalculating normals for 4402 Points, sample size= 88\n",
      "Calculating mesh from points with R = 542.7023788342026\n",
      "Recalculating normals for 4403 Points, sample size= 88\n",
      "Calculating mesh from points with R = 542.6987123320869\n",
      "Recalculating normals for 4391 Points, sample size= 87\n",
      "Calculating mesh from points with R = 542.6965415307668\n",
      "Recalculating normals for 1473 Points, sample size= 29\n",
      "Calculating mesh from points with R = 542.6735932861698\n",
      "Recalculating normals for 4283 Points, sample size= 85\n",
      "Calculating mesh from points with R = 557.7736549330104\n",
      "Recalculating normals for 2455 Points, sample size= 49\n",
      "Calculating mesh from points with R = 348.9759449929405\n",
      "Recalculating normals for 2510 Points, sample size= 50\n",
      "Calculating mesh from points with R = 348.9759449929405\n",
      "Recalculating normals for 2561 Points, sample size= 51\n",
      "Calculating mesh from points with R = 351.4607599877364\n",
      "Recalculating normals for 2272 Points, sample size= 45\n",
      "Calculating mesh from points with R = 348.9759449929405\n",
      "Recalculating normals for 1574 Points, sample size= 31\n",
      "Calculating mesh from points with R = 301.82733317800216\n",
      "Recalculating normals for 2484 Points, sample size= 49\n",
      "Calculating mesh from points with R = 343.9178056904336\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b9e6c9ce676c427f91287ad9c246c868",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Plot(antialias=3, axes=['\\\\sum_i x_i', 'y', 'z'], axes_helper=1.0, background_color=16777215, camera=[6161.490…"
      ]
     },
     "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",
    "\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.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",
    "    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",
    "    points = smoothMLS2D(points, f=0.4)  # smooth cloud, input actor is modified\n",
    "\n",
    "    print(\"Nr of points before cleaning polydata:\", pts1.N())\n",
    "\n",
    "    # impose a min distance among mesh points\n",
    "    points.clean(tol=0.01).legend(\"smooth cloud\")\n",
    "    print(\"             after  cleaning polydata:\", pts1.N())\n",
    "\n",
    "    reco = recoSurface(points, bins=128).legend(\"surf reco\")\n",
    "    plot += reco\n",
    "\n",
    "#filename = 'noddy_out.vtr'\n",
    "#e = load(filename, threshold=50)\n",
    "#value = e.info.keys()\n",
    "#scals = e.scalars('geology')\n",
    "#e.pointColors(scals, cmap='Set3')\n",
    "\n",
    "#plot += e\n",
    "\n",
    "plot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-15., -14., -13., -12., -11., -10.,  -9.,  -8.,  -7.,  -6.,  -5.,\n",
       "        -4.,  -3.,  -2.,  -1.,   0.,   1.,   2.,   3.,   4.,   5.,   6.])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.unique(data['Cat'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0. , 0. , 0.5, 1. ])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "colors[i]"
   ]
  },
  {
   "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
}