2

I have a very large data set (~5 million rows) of X, Y, Z coordinate data in a .dat file. I can import this an plot it as a 3D scatter, but would like to interpolate at 3D surface mesh between the points. I have looked at some examples including this, but can't get it to work without getting an Invalid Index error.

Here's what I have:

from pylab import *
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.mlab import griddata

x, y, z = loadtxt('test.dat', unpack = True)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

xs = x
ys = y
zs = z


fig = plt.figure()
ax = fig.gca(projection='3d')

xi = np.linspace(min(xs), max(xs))
yi = np.linspace(min(ys), max(ys))

X, Y = np.meshgrid(xi, yi)
Z = griddata(xs, ys, zs, xi, yi)

surf = ax.plot_surface(X, Y, Z, rstride=5, cstride=5, cmap=cm.jet,
                       linewidth=1, antialiased=True)

ax.set_zlim3d(np.min(Z), np.max(Z))
fig.colorbar(surf)

plt.show()

And here is the error I get:

python test2.py
Traceback (most recent call last):
  File "test2.py", line 25, in <module>
    Z = griddata(xs, ys, zs, xi, yi)
  File "/Library/Frameworks/Python.framework/Versions/7.3/lib/python2.7/site-packages/matplotlib/mlab.py", line 2768, in griddata
    tri = delaunay.Triangulation(x,y)
  File "/Library/Frameworks/Python.framework/Versions/7.3/lib/python2.7/site-packages/matplotlib/delaunay/triangulate.py", line 90, in __init__
    self.hull = self._compute_convex_hull()
  File "/Library/Frameworks/Python.framework/Versions/7.3/lib/python2.7/site-packages/matplotlib/delaunay/triangulate.py", line 115, in _compute_convex_hull
    edges.update(dict(zip(self.triangle_nodes[border[:,0]][:,1],
IndexError: invalid index

I'm very new to python and can't figure this out. Any ideas? Thanks!

Here is a portion of the data set I am using as a test:

x = [-0.91392505 -0.89382458 -0.87372404 -0.85362357 -0.83352304 -0.81342256
 -0.79332203 -0.77322155 -0.75312102 -0.73302054 -0.71292001 -0.69281954
 -0.672719   -0.65261853 -0.63251805 -0.61241752 -0.59231698 -0.57221651
 -0.55211604 -0.5320155 ]
y = [ 111.25222  111.25222  111.25222  111.25222  111.25222  111.25222
  111.25222  111.25222  111.25222  111.25222  111.25222  111.25222
  111.25222  111.25222  111.25222  111.25222  111.25222  111.25222
  111.25222  111.25222]
z = [  1.42150589e-06   7.77236906e-08   3.33243515e-05   1.70491203e-05
   6.01433214e-08   9.20339880e-06   4.21266122e-06   3.39080143e-04
   1.37568408e-04   7.34613104e-06   1.07246015e-05   3.39363214e-06
   1.09533859e-04   6.55860058e-05   2.83635192e-04   1.87549119e-06
   9.96281233e-06   2.78222607e-04   2.88286283e-05   4.60408737e-05]
Community
  • 1
  • 1
Jonathan Fry
  • 643
  • 2
  • 7
  • 13
  • You might get more answers if you provide some data that demonstrates the problem. – unutbu Mar 28 '13 at 02:01
  • can you re-format your data examples so that they can be copy-pasted? Also, all your Y-values are identical, was that intentional? – tacaswell Apr 01 '13 at 18:33
  • The Y-values being identical is correct. These are only the first 20 values of ~5 million total data points, and the Y values change very slowly over that range. Also, I believe the data are copyable. I just tried and it worked; if not let me know, I'll try and correct it. – Jonathan Fry Apr 01 '13 at 18:57
  • `x` returns a len 1 list which is the sum of everything, and `y` raises an error, you need coma's between the entries. – tacaswell Apr 01 '13 at 19:23

1 Answers1

1

Pretty sure the issue is you are calling griddata wrong (documentation). Try:

Z = scipy.interpolate.griddata((xs, ys), zs, X, Y)
tacaswell
  • 84,579
  • 22
  • 210
  • 199
  • I read up on the documentation for griddata and tried this suggestion. I was not able to get it to plot. I recieved an error saying, "griddata() takes at least 5 arguments (5 given)" or "IndexError: invalid index" – Jonathan Fry Apr 01 '13 at 18:27
  • 1
    `scipy.interpolate.griddata(points, values, xi, method='linear', fill_value=nan)` is not the same as `pylab.griddata(x, y, z, xi, yi, interp='nn')` – szmoore Dec 06 '13 at 04:20
  • @matches Good point. You should post an answer using the mpl version. – tacaswell Dec 06 '13 at 16:04