3d Plots.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context(\"talk\", font_scale=1.5, rc={\"lines.linewidth\": 2.5})\n",
    "# %matplotlib inline\n",
    "\n",
    "# filepath = '/Users/forrest/Documents/CRAWLAB-Student-Code/Forrest Montgomery/CDPM/Motion_Calculations/csv_workspace/'\n",
    "filepath = '/Users/forrest/Documents/CRAWLAB-Student-Code/Forrest Montgomery/CDPM/Motion_Calculations/CSV/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n=20\n",
    "amps = np.array([(0,0)])\n",
    "percent = np.array([0])\n",
    "for i in range(n):\n",
    "    resp = np.genfromtxt(filepath+'response'+str(i)+'.csv', delimiter=',')\n",
    "    shaped = np.genfromtxt(filepath+'response_s'+str(i)+'.csv', delimiter=',')\n",
    "\n",
    "    if shaped[:,0][1] <= 10:\n",
    "        itemindex = np.where(shaped[:,0] >= 10)\n",
    "    else:\n",
    "        itemindex = np.where(shaped[:,0] <= 10)\n",
    "\n",
    "    if resp[:,0][1] <= 10:\n",
    "        itemindex_resp = np.where(resp[:,0] >= 10)\n",
    "    else:\n",
    "        itemindex_resp = np.where(resp[:,0] <= 10)\n",
    "\n",
    "    shaped_ends = itemindex[:1][0][0]\n",
    "\n",
    "    unshaped_ends = itemindex_resp[:1][0][0]\n",
    "\n",
    "    beta_max = np.abs(np.max(resp[:,2][unshaped_ends:]))\n",
    "    beta_min = np.abs(np.min(resp[:,2][unshaped_ends:]))\n",
    "    beta_amp = np.abs(beta_max - beta_min)\n",
    "\n",
    "    beta_shaped_max = np.abs(np.max(shaped[:,2][shaped_ends:]))\n",
    "    beta_shaped_min = np.abs(np.min(shaped[:,2][shaped_ends:]))\n",
    "    beta_shaped_amp = np.abs(beta_shaped_max - beta_shaped_min)\n",
    "\n",
    "    percent_vib = (beta_shaped_amp / beta_amp) * 100\n",
    "    amps = np.append(amps, [(beta_amp, beta_shaped_amp)], axis=0)\n",
    "    percent = np.append(percent, [percent_vib])\n",
    "    percent = np.nan_to_num(percent)\n",
    "amps = amps[1:,:]\n",
    "percent = percent[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def circle(radius, count):\n",
    "    '''\n",
    "    This will take a radius and an amount of points, and create coordinates around\n",
    "    the center point (10, 7.5). So the radius cannot be larger than 7.5.\n",
    "    '''\n",
    "    coord = np.array([(10 + radius, 7.5)])\n",
    "    origin = np.array([(10, 7.5)])\n",
    "    circumfrence = 2 * np.pi * radius\n",
    "    amount = circumfrence / count\n",
    "    n = count\n",
    "    theta = amount / radius\n",
    "    for i in range(n):\n",
    "        y = 7.5 - radius * np.sin(theta*i)\n",
    "        x = 10 + radius * np.cos(theta*i)\n",
    "        coord = np.append([(x,y)], coord, axis=0)\n",
    "    coord = np.reshape(coord, [count+1, 2])\n",
    "    coord = coord[:n]\n",
    "    return coord\n",
    "\n",
    "X_and_Y_points = 20\n",
    "\n",
    "X = np.linspace(1.0, 19.0, num=X_and_Y_points)\n",
    "Y = np.linspace(1.0, 14.0, num=X_and_Y_points)\n",
    "n = X_and_Y_points\n",
    "dims=np.array((0,0))\n",
    "for i in range(n):\n",
    "    x_i = np.repeat(X[i],n)\n",
    "    total = np.column_stack((x_i,Y))\n",
    "    dims = np.vstack((dims,total))\n",
    "dims = dims[1:]\n",
    "\n",
    "############## Variables to Change #########################\n",
    "\n",
    "radius_of_circle = 4\n",
    "\n",
    "number_of_points = 20\n",
    "\n",
    "points = circle(radius_of_circle, number_of_points)\n",
    "# points = dims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "itemindex_resp = np.where(percent >= 100)\n",
    "new = np.delete(percent, itemindex_resp)\n",
    "new_dimsx = np.delete(dims[:,0], itemindex_resp)\n",
    "new_dimsy = np.delete(dims[:,1], itemindex_resp)\n",
    "points = np.column_stack((new_dimsx, new_dimsy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# labels = ['point{0}'.format(i) for i in range(number_of_points)]\n",
    "# plt.scatter(points[:,0], points[:,1])\n",
    "# plt.scatter(10, 7.5)\n",
    "# plt.axes().set_aspect('equal', 'datalim')\n",
    "# for label, x, y in zip(labels, points[:, 0], points[:, 1]):\n",
    "#     plt.annotate(label, xy = (x, y))\n",
    "#         textcoords = 'offset points', ha = 'right', va = 'bottom',\n",
    "#         bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 0.5),\n",
    "#         arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(0)\n",
    "plt.plot(new)\n",
    "# plt.plot(t, shaped[:,0])\n",
    "plt.title('Percent Vibration')\n",
    "plt.xlabel('Move Number')\n",
    "plt.ylabel('Percent Vibration')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(amps[:,0], label='Unshaped')\n",
    "plt.plot(amps[:,1], label='Shaped')\n",
    "plt.title('Total Amplitude after Command')\n",
    "plt.legend()\n",
    "plt.xlabel('Move Number')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "resp = np.genfromtxt(filepath+'response'+str(5)+'.csv', delimiter=',')\n",
    "shaped = np.genfromtxt(filepath+'response_s'+str(5)+'.csv', delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(t, resp[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(t, shaped[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "amps = np.array([(0,0)])\n",
    "percent = np.array([0])\n",
    "if shaped[:,0][1] >= 10:\n",
    "    itemindex = np.where(shaped[:,0] <= 10)\n",
    "else:\n",
    "    itemindex = np.where(shaped[:,0] >= 10)\n",
    "\n",
    "if resp[:,0][1] >= 10:\n",
    "    itemindex_resp = np.where(resp[:,0] <= 10)\n",
    "else:\n",
    "    itemindex_resp = np.where(resp[:,0] >= 10)\n",
    "\n",
    "shaped_ends = itemindex[:1][0][0]\n",
    "\n",
    "unshaped_ends = itemindex_resp[:1][0][0]\n",
    "\n",
    "beta_max = np.abs(np.max(resp[:,2][unshaped_ends:]))\n",
    "beta_min = np.abs(np.min(resp[:,2][unshaped_ends:]))\n",
    "beta_amp = np.abs(beta_max - beta_min)\n",
    "\n",
    "beta_shaped_max = np.abs(np.max(shaped[:,2][shaped_ends:]))\n",
    "beta_shaped_min = np.abs(np.min(shaped[:,2][shaped_ends:]))\n",
    "beta_shaped_amp = np.abs(beta_shaped_max - beta_shaped_min)\n",
    "\n",
    "percent_vib = (beta_shaped_amp / beta_amp) * 100\n",
    "amps = np.append(amps, [(beta_amp, beta_shaped_amp)], axis=0)\n",
    "percent = np.append(percent, [percent_vib])\n",
    "percent = np.nan_to_num(percent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(shaped[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "percent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# import numpy as np\n",
    "# from matplotlib.pyplot import *             # Grab MATLAB plotting functions\n",
    "# from mpl_toolkits.mplot3d import Axes3D     # 3D plotting functions\n",
    "# from scipy.interpolate import griddata      # interpolate package's griddate function is better\n",
    "\n",
    "\n",
    "# # Example data so that code will show a plot... \n",
    "# # It's best to remove and fill in your own\n",
    "# x = points[:,0]\n",
    "# y = points[:,1]\n",
    "\n",
    "# z = new\n",
    "    \n",
    "\n",
    "# # This function will take x and y arrays and create matrices from them\n",
    "# #   If X is shape (N,) or (N,1) and y is shape (M,) or (M,1), then \n",
    "# #   x_grid and y_grid will be (M, N)\n",
    "# x_grid, y_grid = np.meshgrid(x, y)\n",
    "\n",
    "\n",
    "# heights = griddata((x, y),z,(x_grid, y_grid),method='nearest')\n",
    "\n",
    "# fig = figure(figsize=(9,6))\n",
    "# subplots_adjust(bottom=0.17,left=0.17,top=0.96,right=0.96)\n",
    "# ax = Axes3D(fig)\n",
    "\n",
    "\n",
    "# surf = ax.plot_surface(x_grid, y_grid, heights, linewidth=0, alpha=0.85,cmap=cm.bwr)\n",
    "# # cset = ax.contourf(x_grid, y_grid, heights, zdir='z', offset= 0, cmap=cm.bwr)\n",
    "# # cset = ax.contour(x_grid, y_grid, heights, zdir='x', offset= np.max(x), cmap=cm.gray)\n",
    "# # cset = ax.contour(x_grid, y_grid, heights, zdir='y', offset=1.5*np.max(y), cmap=cm.gray)\n",
    "\n",
    "\n",
    "# # adjusts the padding around the 3D plot\n",
    "# ax.dist = 11\n",
    "\n",
    "# # Change the colorbar font\n",
    "# color_bar = colorbar(surf,shrink=0.5,aspect=8)\n",
    "# cbytick_obj = getp(color_bar.ax.axes, 'yticklabels')\n",
    "# setp(cbytick_obj, family='CMUSerif-Roman',fontsize=18)\n",
    "\n",
    "# # Export to svg for final edits in iDraw\n",
    "# # savefig('3D_plot.svg')\n",
    "# show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# data = np.genfromtxt('300.txt')\n",
    "x = points[:,0]\n",
    "y = points[:,1]\n",
    "z = new\n",
    "p = ax.scatter(x, y, z, c='r', marker='o')\n",
    "\n",
    "ax.set_xlabel('X Label')\n",
    "ax.set_ylabel('Y Label')\n",
    "ax.set_zlabel('Z Label')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Argument 'zs' must be of same size as 'xs' and 'ys' or of size 1.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-9003e07e3b47>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;31m# for i in range(20):\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;31m# ax.plot(0, 0, 10)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0msc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvmin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvmax\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m35\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_xlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'X (m)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_ylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Y (m)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m//anaconda/lib/python3.5/site-packages/mpl_toolkits/mplot3d/axes3d.py\u001b[0m in \u001b[0;36mscatter\u001b[0;34m(self, xs, ys, zs, zdir, s, c, depthshade, *args, **kwargs)\u001b[0m\n\u001b[1;32m   2258\u001b[0m                 \u001b[0mzs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mzs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2259\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2260\u001b[0;31m                 raise ValueError((\"Argument 'zs' must be of same size as 'xs' \"\n\u001b[0m\u001b[1;32m   2261\u001b[0m                     \"and 'ys' or of size 1.\"))\n\u001b[1;32m   2262\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: Argument 'zs' must be of same size as 'xs' and 'ys' or of size 1."
     ]
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "plt.gca().patch.set_facecolor('white')\n",
    "import matplotlib.pyplot as plt\n",
    "cm = plt.cm.get_cmap('inferno')\n",
    "x = points[:,0]\n",
    "y = points[:,1]\n",
    "z = new\n",
    "# for i in range(20):\n",
    "# ax.plot(0, 0, 10)\n",
    "sc = ax.scatter(x, y, z, c=z, vmin=0, vmax=100, s=35, cmap=cm)\n",
    "ax.set_xlabel('X (m)')\n",
    "ax.set_ylabel('Y (m)')\n",
    "ax.set_zlabel('Percent Vibration')\n",
    "ax.set_zlim(0, 35)\n",
    "plt.colorbar(sc)\n",
    "# ax.view_init(elev=9, azim=-137)\n",
    "# plt.savefig('Files.svg', format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# fig = plt.figure(figsize=(6,6))\n",
    "# ax = fig.add_subplot(111)\n",
    "# ax.set_title(\"X vs Y AVG\",fontsize=14)\n",
    "# ax.set_xlabel(\"XAVG\",fontsize=12)\n",
    "# ax.set_ylabel(\"YAVG\",fontsize=12)\n",
    "# ax.grid(True,linestyle='-',color='0.75')\n",
    "# x = points[:,0]\n",
    "# y = points[:,1]\n",
    "# z = new\n",
    "# plt.axes().set_aspect('equal', 'datalim')\n",
    "# # scatter with colormap mapping to z value\n",
    "# ax.scatter(x,y,s=20,c=z, marker = 'o', cmap = cm.jet );\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "x, y and z must be equal-length 1-D arrays",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-bfac1f603aaa>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmeshgrid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0myi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m \u001b[0mZ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgriddata\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0myi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap=cm.jet,\n",
      "\u001b[0;32m//anaconda/lib/python3.5/site-packages/matplotlib/mlab.py\u001b[0m in \u001b[0;36mgriddata\u001b[0;34m(x, y, z, xi, yi, interp)\u001b[0m\n\u001b[1;32m   3385\u001b[0m     \u001b[0mz\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masanyarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3386\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3387\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"x, y and z must be equal-length 1-D arrays\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3388\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3389\u001b[0m     \u001b[0mxi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masanyarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: x, y and z must be equal-length 1-D arrays"
     ]
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.mlab import griddata\n",
    "import numpy as np\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "# data = np.genfromtxt('300.txt')\n",
    "x = points[:,0]\n",
    "y = points[:,1]\n",
    "z = new\n",
    "\n",
    "xi = np.linspace(min(x), max(x))\n",
    "yi = np.linspace(min(y), max(y))\n",
    "\n",
    "X, Y = np.meshgrid(xi, yi)\n",
    "Z = griddata(x, y, z, xi, yi)\n",
    "\n",
    "surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap=cm.jet,\n",
    "                       linewidth=1, antialiased=True)\n",
    "\n",
    "surf.set_clim(vmin=0, vmax=100)\n",
    "ax.set_zlim3d(0, 100)\n",
    "ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "fig.colorbar(surf)\n",
    "ax.text2D(0.05, 0.95, \"2D Text\", transform=ax.transAxes)\n",
    "ax.set_xlabel('X axis')\n",
    "ax.set_ylabel('Y axis')\n",
    "ax.set_zlabel('Z axis')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import numpy as np\n",
    "# import matplotlib.pyplot as plt\n",
    "# import matplotlib.mlab as ml\n",
    "\n",
    "# # np.random.seed(8)\n",
    "# ndata = 10\n",
    "# # ny, nx = 100, 200\n",
    "\n",
    "# x = points[:,0]\n",
    "# y = points[:,1]\n",
    "# z = percent\n",
    "# xmin, xmax = np.min(x), np.max(x)\n",
    "# ymin, ymax = np.min(y), np.max(y)\n",
    "\n",
    "# # x = np.linspace(1, 10, ndata)\n",
    "# # y = np.linspace(1, 10, ndata)\n",
    "# # z = np.random.random(ndata)\n",
    "# x = np.r_[x,xmin,xmax]\n",
    "# y = np.r_[y,ymax,ymin]\n",
    "# z = np.r_[z,z[0],z[-1]]\n",
    "# xi = np.linspace(xmin, xmax, nx)\n",
    "# yi = np.linspace(ymin, ymax, ny)\n",
    "\n",
    "\n",
    "# # Requires installation of natgrid\n",
    "# # http://sourceforge.net/projects/matplotlib/files/matplotlib-toolkits/\n",
    "# zi = ml.griddata(x, y, z, xi, yi, interp='nn')\n",
    "\n",
    "# # Or, without natgrid:\n",
    "# # zi = ml.griddata(x, y, z, xi, yi, interp='linear')\n",
    "\n",
    "# # plt.contour(xi, yi, zi, 15, linewidths = 0.5, colors = 'k')\n",
    "# plt.pcolormesh(xi, yi, zi, cmap = plt.get_cmap('rainbow'))\n",
    "\n",
    "# plt.colorbar() \n",
    "# plt.scatter(x, y, marker = 'o', c = 'b', s = 10, zorder = 10)\n",
    "# # plt.xlim(xmin, xmax)\n",
    "# # plt.ylim(ymin, ymax)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# from mpl_toolkits.mplot3d import Axes3D\n",
    "# from matplotlib import cm\n",
    "# import matplotlib.pyplot as plt\n",
    "# import numpy as np\n",
    "\n",
    "# n_angles = 36\n",
    "# n_radii = 8\n",
    "\n",
    "# # An array of radii\n",
    "# # Does not include radius r=0, this is to eliminate duplicate points\n",
    "# radii = np.linspace(0.125, 1.0, n_radii)\n",
    "\n",
    "# # An array of angles\n",
    "# angles = np.linspace(0, 2*np.pi, n_angles, endpoint=False)\n",
    "\n",
    "# # Repeat all angles for each radius\n",
    "# angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)\n",
    "\n",
    "# # Convert polar (radii, angles) coords to cartesian (x, y) coords\n",
    "# # (0, 0) is added here. There are no duplicate points in the (x, y) plane\n",
    "# x = np.append(0, (radii*np.cos(angles)).flatten())\n",
    "# y = np.append(0, (radii*np.sin(angles)).flatten())\n",
    "\n",
    "# # Pringle surface\n",
    "# z = np.sin(-x*y)\n",
    "\n",
    "# fig = plt.figure()\n",
    "# ax = fig.gca(projection='3d')\n",
    "\n",
    "# ax.plot_trisurf(x,y,z, cmap=cm.jet, linewidth=0.2)\n",
    "\n",
    "# plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.mlab import griddata\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "import numpy as np\n",
    "from scipy.stats import gaussian_kde\n",
    "from matplotlib import rc\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "rc('font',size=28)\n",
    "rc('font',family='serif')\n",
    "rc('axes',labelsize=32)\n",
    "\n",
    "# filename = sys.argv[1]\n",
    "# data = np.load(filename)\n",
    "# output = data[171,:,:].flatten()*np.sqrt(2.0)/13074\n",
    "x = points[:,0]\n",
    "y = points[:,1]\n",
    "z = new\n",
    "\n",
    "bins=30\n",
    "fig = plt.figure(figsize=(13.5,8))  # PRL default width\n",
    "\n",
    "ax1 = fig.add_subplot(121, projection='3d')\n",
    "plt.subplots_adjust(left=0, right=0.90, top=1, bottom=0, wspace=0.22)\n",
    "\n",
    "hist, xedges, yedges = np.histogram2d(x, y, bins)\n",
    "X, Y = np.meshgrid(xedges[:-1], yedges[:-1])\n",
    "Z = hist\n",
    "Z = Z/Z.sum() # normalize Z\n",
    "surf = ax1.plot_surface(X,Y,Z,cstride=1,rstride=1,color=\"white\",shade=False)\n",
    "contourz = (Z.min()-Z.max())*0.4  # where to put the contours\n",
    "cset = ax1.contour(X,Y,Z.reshape(X.shape),zdir='z',offset=contourz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ax2 = fig.add_subplot(122, projection='3d')\n",
    "hist, xedges, yedges = np.histogram2d(x, y, bins)\n",
    "X, Y = np.meshgrid(xedges[:-1], yedges[:-1])\n",
    "Z = hist\n",
    "Z = Z/Z.sum() # normalize Z\n",
    "surf2 = ax2.plot_surface(X,Y,Z,cstride=1,rstride=1,color=\"white\",shade=False)\n",
    "contourz2 = (Z.min()-Z.max())*0.4  # where to put the contours\n",
    "cset2 = ax2.contour(X,Y,Z.reshape(X.shape),zdir='z',offset=contourz2)\n",
    "\n",
    "ax1.set_xlabel(r'$x_p$')\n",
    "ax1.set_ylabel(r'$y_p$')\n",
    "#ax1.set_zlabel(r'$Q(x_p,y_p)$')\n",
    "ax1.text(X.min()*1.1, Y.min()*1.1, Z.max()*1.1, r'$Q(x_p,y_p)$')\n",
    "ax1.xaxis.set_rotate_label(False)\n",
    "ax1.yaxis.set_rotate_label(False)\n",
    "ax1.zaxis.set_rotate_label(False)\n",
    "ax1.set_zlim(contourz,Z.max()*1.1)\n",
    "ax1.view_init(elev=10, azim=135)\n",
    "ax1.grid(False)\n",
    "ax1.xaxis.pane.set_edgecolor('black')\n",
    "ax1.yaxis.pane.set_edgecolor('black')\n",
    "ax1.xaxis.pane.fill = False\n",
    "ax1.yaxis.pane.fill = False\n",
    "ax1.zaxis.pane.fill = False\n",
    "ax1.xaxis.set_major_locator(MultipleLocator(5))\n",
    "ax1.yaxis.set_major_locator(MultipleLocator(5))\n",
    "ax1.zaxis.set_major_locator(MultipleLocator(0.01))\n",
    "[t.set_va('center') for t in ax1.get_yticklabels()]\n",
    "[t.set_ha('left') for t in ax1.get_yticklabels()]\n",
    "[t.set_va('center') for t in ax1.get_xticklabels()]\n",
    "[t.set_ha('right') for t in ax1.get_xticklabels()]\n",
    "[t.set_va('center') for t in ax1.get_zticklabels()]\n",
    "[t.set_ha('left') for t in ax1.get_zticklabels()]\n",
    "ax1.xaxis._axinfo['tick']['inward_factor'] = 0\n",
    "ax1.xaxis._axinfo['tick']['outward_factor'] = 0.4\n",
    "ax1.yaxis._axinfo['tick']['inward_factor'] = 0\n",
    "ax1.yaxis._axinfo['tick']['outward_factor'] = 0.4\n",
    "ax1.zaxis._axinfo['tick']['inward_factor'] = 0\n",
    "ax1.zaxis._axinfo['tick']['outward_factor'] = 0.4\n",
    "ax1.zaxis._axinfo['tick']['outward_factor'] = 0.4\n",
    "\n",
    "ax2.set_xlabel(r'$x_p$',fontsize=32)\n",
    "ax2.set_ylabel(r'$y_p$')\n",
    "#ax2.set_zlabel(r'$Q(x_p,y_p)$')\n",
    "ax2.text(X.min()*1.1, Y.min()*1.1, Z.max()*1.1, r'$Q(x_p,y_p)$')\n",
    "ax2.xaxis.set_rotate_label(False)\n",
    "ax2.yaxis.set_rotate_label(False)\n",
    "ax2.zaxis.set_rotate_label(False)\n",
    "ax2.set_zlim(contourz,Z.max()*1.1)\n",
    "ax2.view_init(elev=10, azim=135)\n",
    "ax2.grid(False)\n",
    "ax2.xaxis.pane.set_edgecolor('black')\n",
    "ax2.yaxis.pane.set_edgecolor('black')\n",
    "ax2.xaxis.pane.fill = False\n",
    "ax2.yaxis.pane.fill = False\n",
    "ax2.zaxis.pane.fill = False\n",
    "ax2.xaxis.set_major_locator(MultipleLocator(5))\n",
    "ax2.yaxis.set_major_locator(MultipleLocator(5))\n",
    "ax2.zaxis.set_major_locator(MultipleLocator(0.01))\n",
    "[t.set_va('center') for t in ax2.get_yticklabels()]\n",
    "[t.set_ha('left') for t in ax2.get_yticklabels()]\n",
    "[t.set_va('center') for t in ax2.get_xticklabels()]\n",
    "[t.set_ha('right') for t in ax2.get_xticklabels()]\n",
    "[t.set_va('center') for t in ax2.get_zticklabels()]\n",
    "[t.set_ha('left') for t in ax2.get_zticklabels()]\n",
    "ax2.xaxis._axinfo['tick']['inward_factor'] = 0\n",
    "ax2.xaxis._axinfo['tick']['outward_factor'] = 0.4\n",
    "ax2.yaxis._axinfo['tick']['inward_factor'] = 0\n",
    "ax2.yaxis._axinfo['tick']['outward_factor'] = 0.4\n",
    "ax2.zaxis._axinfo['tick']['inward_factor'] = 0\n",
    "ax2.zaxis._axinfo['tick']['outward_factor'] = 0.4\n",
    "\n",
    "#ax.set_xlim(-40, 40)\n",
    "#ax.set_ylim(-40, 40)\n",
    "#ax.zaxis.set_major_locator(LinearLocator(3))\n",
    "#ax.zaxis.set_major_formatter(FormatStrFormatter(\"%.02f\"))\n",
    "#ax.locator_params(tight=True)\n",
    "#ax.tick_params(labelsize=30)\n",
    "#ax.set_xticks([min(x), max(x), 0.0])\n",
    "#ax.set_yticks([min(y), max(y), 0.0])\n",
    "#ax.set_zticks([min(Z), max(Z), 0.0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "different number of values and points",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-b00e28baefe1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m \u001b[0;31m# Interpolate; there's also method='cubic' for 2-D data such as here\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 14\u001b[0;31m \u001b[0mzi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgriddata\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mxi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0myi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'linear'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0mcm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_cmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'RdYlBu'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m//anaconda/lib/python3.5/site-packages/scipy/interpolate/ndgriddata.py\u001b[0m in \u001b[0;36mgriddata\u001b[0;34m(points, values, xi, method, fill_value, rescale)\u001b[0m\n\u001b[1;32m    215\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'linear'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    216\u001b[0m         ip = LinearNDInterpolator(points, values, fill_value=fill_value,\n\u001b[0;32m--> 217\u001b[0;31m                                   rescale=rescale)\n\u001b[0m\u001b[1;32m    218\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    219\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'cubic'\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mndim\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mscipy/interpolate/interpnd.pyx\u001b[0m in \u001b[0;36mscipy.interpolate.interpnd.LinearNDInterpolator.__init__ (scipy/interpolate/interpnd.c:4934)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mscipy/interpolate/interpnd.pyx\u001b[0m in \u001b[0;36mscipy.interpolate.interpnd.NDInterpolatorBase.__init__ (scipy/interpolate/interpnd.c:2386)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mscipy/interpolate/interpnd.pyx\u001b[0m in \u001b[0;36mscipy.interpolate.interpnd._check_init_shape (scipy/interpolate/interpnd.c:4556)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: different number of values and points"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.interpolate\n",
    "\n",
    "# Generate data: for N=1e6, the triangulation hogs 1 GB of memory\n",
    "x = points[:,0]\n",
    "y = points[:,1]\n",
    "z = new\n",
    "# Set up a regular grid of interpolation points\n",
    "xi, yi = np.linspace(x.min(), x.max(), 300), np.linspace(y.min(), y.max(), 300)\n",
    "xi, yi = np.meshgrid(xi, yi)\n",
    "\n",
    "# Interpolate; there's also method='cubic' for 2-D data such as here\n",
    "zi = scipy.interpolate.griddata((x, y), z, (xi, yi), method='linear')\n",
    "\n",
    "cm = plt.cm.get_cmap('RdYlBu')\n",
    "plt.imshow(zi, vmin=z.min(), vmax=z.max(), origin='lower',cmap=cm,\n",
    "           extent=[x.min(), x.max(), y.min(), y.max()])\n",
    "plt.colorbar()\n",
    "plt.ylim(15,0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def griddata(x, y, z, binsize=0.01, retbin=True, retloc=True):\n",
    "    \"\"\"\n",
    "    Place unevenly spaced 2D data on a grid by 2D binning (nearest\n",
    "    neighbor interpolation).\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : ndarray (1D)\n",
    "        The idependent data x-axis of the grid.\n",
    "    y : ndarray (1D)\n",
    "        The idependent data y-axis of the grid.\n",
    "    z : ndarray (1D)\n",
    "        The dependent data in the form z = f(x,y).\n",
    "    binsize : scalar, optional\n",
    "        The full width and height of each bin on the grid.  If each\n",
    "        bin is a cube, then this is the x and y dimension.  This is\n",
    "        the step in both directions, x and y. Defaults to 0.01.\n",
    "    retbin : boolean, optional\n",
    "        Function returns `bins` variable (see below for description)\n",
    "        if set to True.  Defaults to True.\n",
    "    retloc : boolean, optional\n",
    "        Function returns `wherebins` variable (see below for description)\n",
    "        if set to True.  Defaults to True.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    grid : ndarray (2D)\n",
    "        The evenly gridded data.  The value of each cell is the median\n",
    "        value of the contents of the bin.\n",
    "    bins : ndarray (2D)\n",
    "        A grid the same shape as `grid`, except the value of each cell\n",
    "        is the number of points in that bin.  Returns only if\n",
    "        `retbin` is set to True.\n",
    "    wherebin : list (2D)\n",
    "        A 2D list the same shape as `grid` and `bins` where each cell\n",
    "        contains the indicies of `z` which contain the values stored\n",
    "        in the particular bin.\n",
    "\n",
    "    Revisions\n",
    "    ---------\n",
    "    2010-07-11  ccampo  Initial version\n",
    "    \"\"\"\n",
    "    # get extrema values.\n",
    "    xmin, xmax = x.min(), x.max()\n",
    "    ymin, ymax = y.min(), y.max()\n",
    "\n",
    "    # make coordinate arrays.\n",
    "    xi      = np.arange(xmin, xmax+binsize, binsize)\n",
    "    yi      = np.arange(ymin, ymax+binsize, binsize)\n",
    "    xi, yi = np.meshgrid(xi,yi)\n",
    "\n",
    "    # make the grid.\n",
    "    grid           = np.zeros(xi.shape, dtype=x.dtype)\n",
    "    nrow, ncol = grid.shape\n",
    "    if retbin: bins = np.copy(grid)\n",
    "\n",
    "    # create list in same shape as grid to store indices\n",
    "    if retloc:\n",
    "        wherebin = np.copy(grid)\n",
    "        wherebin = wherebin.tolist()\n",
    "\n",
    "    # fill in the grid.\n",
    "    for row in range(nrow):\n",
    "        for col in range(ncol):\n",
    "            xc = xi[row, col]    # x coordinate.\n",
    "            yc = yi[row, col]    # y coordinate.\n",
    "\n",
    "            # find the position that xc and yc correspond to.\n",
    "            posx = np.abs(x - xc)\n",
    "            posy = np.abs(y - yc)\n",
    "            ibin = np.logical_and(posx < binsize/2., posy < binsize/2.)\n",
    "            ind  = np.where(ibin == True)[0]\n",
    "\n",
    "            # fill the bin.\n",
    "            bin = z[ibin]\n",
    "            if retloc: wherebin[row][col] = ind\n",
    "            if retbin: bins[row, col] = bin.size\n",
    "            if bin.size != 0:\n",
    "                binval         = np.median(bin)\n",
    "                grid[row, col] = binval\n",
    "            else:\n",
    "                grid[row, col] = np.nan   # fill empty bins with nans.\n",
    "\n",
    "    # return the grid\n",
    "    if retbin:\n",
    "        if retloc:\n",
    "            return grid, bins, wherebin\n",
    "        else:\n",
    "            return grid, bins\n",
    "    else:\n",
    "        if retloc:\n",
    "            return grid, wherebin\n",
    "        else:\n",
    "            return grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib.mlab import griddata\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#-- Read the data.\n",
    "# I'm going to use `pandas` to read in and work with your data, mostly due to\n",
    "# the text site names. Using pandas is optional, however.\n",
    "x = points[:,0]\n",
    "y = points[:,1]\n",
    "z = new\n",
    "\n",
    "#-- Now let's grid your data.\n",
    "# First we'll make a regular grid to interpolate onto. This is equivalent to\n",
    "# your call to `mgrid`, but it's broken down a bit to make it easier to\n",
    "# understand. The \"30j\" in mgrid refers to 30 rows or columns.\n",
    "numcols, numrows = 30, 30\n",
    "xi = np.linspace(x.min(), x.max(), numcols)\n",
    "yi = np.linspace(y.min(), y.max(), numrows)\n",
    "xi, yi = np.meshgrid(xi, yi)\n",
    "\n",
    "#-- Interpolate at the points in xi, yi\n",
    "# \"griddata\" expects \"raw\" numpy arrays, so we'll pass in\n",
    "# data.x.values instead of just the pandas series data.x\n",
    "zi = griddata(x, y, z, xi, yi)\n",
    "\n",
    "#-- Display the results\n",
    "fig, ax = plt.subplots()\n",
    "im = ax.contourf(xi, yi, zi)\n",
    "ax.scatter(x,y, c=z, s=100,\n",
    "           vmin=zi.min(), vmax=zi.max())\n",
    "fig.colorbar(im)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}