title: What to do if I want 3D spline/smooth interpolation of random unstructured data?
body: <p>I was inspired by this <a href="https://stackoverflow.com/a/28506619/3904031">answer</a> by @James to see how <code>griddata</code> and <code>map_coordinates</code> might be used. In the examples below I'm showing 2D data, but my interest is in 3D. I noticed that <code>griddata</code> only provides splines for 1D and 2D, and is limited to linear interpolation for 3D and higher (probably for very good reasons). However, map_coordinates seems to be fine with 3D using higher order (smoother than piece-wise linear) interpolation.</p>
  
  <p><strong>My primary question:</strong>  if I have random, unstructured data (where I can not use map_coordinates) in 3D, is there some way to get smoother than piece-wise linear interpolation within the NumPy SciPy universe, or at least nearby? </p>
  
  <p>My secondary question: is spline for 3D not available in <code>griddata</code> because it is difficult or tedious to implement, or is there a fundamental difficulty?</p>
  
  <p>The images and horrible python below show my current understanding of how griddata and map_coordinates can or can't be used. Interpolation is done along the thick black line.</p>
  
  <p><strong>STRUCTURED DATA:</strong></p>
  
  <p><a href="https://i.stack.imgur.com/fB9nZ.png" rel="noreferrer"><img src="https://i.stack.imgur.com/fB9nZ.png" alt="interpolation of structured data"></a></p>
  
  <p><strong><em>UN</em>STRUCTURED DATA:</strong></p>
  
  <p><a href="https://i.stack.imgur.com/od6x4.png" rel="noreferrer"><img src="https://i.stack.imgur.com/od6x4.png" alt="interpolation of unstructured data"></a></p>
  
  <p><strong>Horrible python:</strong></p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  
  def g(x, y):
      return np.exp(-((x-1.0)**2 + (y-1.0)**2))
  
  def findit(x, X):  # or could use some 1D interpolation
      fraction = (x - X[0]) / (X[-1]-X[0])
      return fraction * float(X.shape[0]-1)
  
  nth, nr = 12, 11
  theta_min, theta_max = 0.2, 1.3
  r_min,     r_max     = 0.7, 2.0
  
  theta = np.linspace(theta_min, theta_max, nth)
  r     = np.linspace(r_min,     r_max,     nr)
  
  R, TH = np.meshgrid(r, theta)
  Xp, Yp  = R*np.cos(TH), R*np.sin(TH)
  array = g(Xp, Yp)
  
  x, y = np.linspace(0.0, 2.0, 200), np.linspace(0.0, 2.0, 200)
  X, Y = np.meshgrid(x, y)
  blob = g(X, Y)
  
  xtest = np.linspace(0.25, 1.75, 40)
  ytest = np.zeros_like(xtest) + 0.75
  
  rtest = np.sqrt(xtest**2 + ytest**2)
  thetatest = np.arctan2(xtest, ytest)
  
  ir = findit(rtest, r)
  it = findit(thetatest, theta)
  
  plt.figure()
  
  plt.subplot(2,1,1)
  
  plt.scatter(100.0*Xp.flatten(), 100.0*Yp.flatten())
  plt.plot(100.0*xtest, 100.0*ytest, '-k', linewidth=3)
  plt.hold
  plt.imshow(blob, origin='lower', cmap='gray')
  
  plt.text(5, 5, "don't use jet!", color='white')
  
  
  exact = g(xtest, ytest)
  
  import scipy.ndimage.interpolation as spndint
  ndint0 = spndint.map_coordinates(array, [it, ir], order=0)
  ndint1 = spndint.map_coordinates(array, [it, ir], order=1)
  ndint2 = spndint.map_coordinates(array, [it, ir], order=2)
  
  import scipy.interpolate as spint
  points = np.vstack((Xp.flatten(), Yp.flatten())).T   # could use np.array(zip(...))
  grid_x = xtest
  grid_y = np.array([0.75])
  
  g0 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='nearest')
  g1 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='linear')
  g2 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='cubic')
  
  
  plt.subplot(4,2,5)
  
  plt.plot(exact, 'or')
  #plt.plot(ndint0)
  plt.plot(ndint1)
  plt.plot(ndint2)
  plt.title("map_coordinates")
  
  plt.subplot(4,2,6)
  
  plt.plot(exact, 'or')
  #plt.plot(g0)
  plt.plot(g1)
  plt.plot(g2)
  plt.title("griddata")
  
  plt.subplot(4,2,7)
  
  #plt.plot(ndint0 - exact)
  plt.plot(ndint1 - exact)
  plt.plot(ndint2 - exact)
  plt.title("error map_coordinates")
  
  plt.subplot(4,2,8)
  
  #plt.plot(g0 - exact)
  plt.plot(g1 - exact)
  plt.plot(g2 - exact)
  plt.title("error griddata")
  
  plt.show()
  
  
  seed_points_rand = 2.0 * np.random.random((400, 2))
  rr = np.sqrt((seed_points_rand**2).sum(axis=-1))
  thth = np.arctan2(seed_points_rand[...,1], seed_points_rand[...,0])
  isinside = (rr&gt;r_min) * (rr&lt;r_max) * (thth&gt;theta_min) * (thth&lt;theta_max)
  points_rand = seed_points_rand[isinside]
  
  Xprand, Yprand = points_rand.T  # unpack
  array_rand = g(Xprand, Yprand)
  
  grid_x = xtest
  grid_y = np.array([0.75])
  
  plt.figure()
  
  plt.subplot(2,1,1)
  
  plt.scatter(100.0*Xprand.flatten(), 100.0*Yprand.flatten())
  plt.plot(100.0*xtest, 100.0*ytest, '-k', linewidth=3)
  plt.hold
  plt.imshow(blob, origin='lower', cmap='gray')
  plt.text(5, 5, "don't use jet!", color='white')
  
  
  g0rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='nearest')
  g1rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='linear')
  g2rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='cubic')
  
  plt.subplot(4,2,6)
  
  plt.plot(exact, 'or')
  #plt.plot(g0rand)
  plt.plot(g1rand)
  plt.plot(g2rand)
  plt.title("griddata")
  
  
  plt.subplot(4,2,8)
  
  #plt.plot(g0rand - exact)
  plt.plot(g1rand - exact)
  plt.plot(g2rand - exact)
  plt.title("error griddata")
  
  plt.show()
  </code></pre>