2

I'm trying to convert MATLAB code to equivalent python. I have 3 arrays and I want to compute interp2d:

nuA = np.asarray([2.439,2.5,2.6,2.7,2.8,3.0,3.2,3.5,4.0,5.0,6.0,8.0,10,15,25])
nuB = np.asarray([0,0.1,0.2,0.3,0.5,0.7,1])
a, b = np.meshgrid(nuA, nuB)
betaTab  = np.transpose(np.asarray([[0.0,2.16,1.0,1.0,1.0,1.0,1.0],[0.0,1.592,3.39,1.0,1.0,1.0,1.0],[0.0,0.759,1.8,1.0,1.0,1.0,1.0],[0.0,0.482,1.048,1.694,1.0,1.0,1.0],[0.0,0.36,0.76,1.232,2.229,1.0,1.0],[0.0,0.253,0.518,0.823,1.575,1.0,1.0],[0.0,0.203,0.41,0.632,1.244,1.906,1.0],[0.0,0.165,0.332,0.499,0.943,1.56,1.0],[0.0,0.136,0.271,0.404,0.689,1.23,2.195],[0.0,0.109,0.216,0.323,0.539,0.827,1.917],[0.0,0.096,0.19,0.284,0.472,0.693,1.759],[0.0,0.082,0.163,0.243,0.412,0.601,1.596],[0.0,0.074,0.147,0.22,0.377,0.546,1.482],[0.0,0.064,0.128,0.191,0.33,0.478,1.362],[0.0,0.056,0.112,0.167,0.285,0.428,1.274]]))
ip = scipy.interpolate.interp2d(a,b,betaTab)

when I try to run it, this warning is displayed:

/usr/local/lib/python2.7/dist-packages/scipy/interpolate/fitpack.py:981: RuntimeWarning: No more knots can be added because the additional knot would
coincide with an old one. Probable cause: s too small or too large
a weight to an inaccurate data point. (fp>s)
    kx,ky=1,1 nx,ny=4,14 m=105 fp=21.576347 s=0.000000
  warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess))

I know that interp2d is different from matlab interp2 and in python RectBivariateSpline function is preferred. But I can't use the latter function because of the length of my data. Also, the final result of ip(xi,yi) is different from the MATLAB answer.

How can I compute interp2d without warning and compute it correctly?

shozdeh
  • 172
  • 1
  • 11
  • What kind of interpolation is the MATLAB code doing? splines, piecewise linear (triangular planes)? – hpaulj Jan 15 '16 at 22:08
  • The error is, in effect, saying that for your data, a linear spline interpolation is bad. It might, for example, be too flat. What does the surface look like? – hpaulj Jan 15 '16 at 22:17

1 Answers1

3

Your input data seems to be quite ill-defined. Here's a surface plot of your input points:

plot of input data, quite ugly

This is not an easy problem to interpolate. Incidentally, I've recently ran into problems where interp2d couldn't even interpolate a smooth data set. So I would suggest checking out scipy.interpolate.griddata instead:

import numpy as np
import scipy.interpolate as interp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#define your data as you did in your question: a, b and betaTab

ip = interp.interp2d(a,b,betaTab)           # original interpolator

aplotv = np.linspace(a.min(),a.max(),100)   # to interpolate at
bplotv = np.linspace(b.min(),b.max(),100)   # to interpolate at
aplot,bplot = np.meshgrid(aplotv,bplotv)    # mesh to interpolate at

# actual values from interp2d:
betainterp2d = ip(aplotv,bplotv)

# actual values from griddata:
betagriddata = interp.griddata(np.array([a.ravel(),b.ravel()]).T,betaTab.ravel(),np.array([aplot.ravel(),bplot.ravel()]).T)
# ^ this probably could be written in a less messy way,
# I'll keep thinking about it

#plot results
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(aplot,bplot,betainterp2d,cmap='viridis',cstride=1,rstride=1)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(aplot,bplot,betagriddata,cmap='viridis',cstride=1,rstride=1)

Results: (left: interp2d, right: griddata)

result interp2d'd result griddata'd

Conclusion: use scipy.interpolate.griddata.

Community
  • 1
  • 1
  • How can I get value from specific point? Result for `griddata` is a 1D `(shape=(10000,))` array while in `ip = interp2d(...)` we can get final interpolation value using `ip(xi,yi)`. Corresponding values of `(xi , yi)` in one scenario are `(2.439, 0.097)` and result must be about `2.1`. – shozdeh Jan 16 '16 at 07:16
  • Solved: Using `print(griddata(np.array([a.ravel(),b.ravel()]).T,betaTab.ravel(),np.array([2.439, 0.097]).T))` get `2.09`. – shozdeh Jan 16 '16 at 07:19
  • That is correct. While an interpolator function from `interp2d` accepts points on a 1d mesh (and outputs interpolated points on a 2d mesh generated from those 1d meshes), `griddata` goes point-wise over the input (and both input and output points should be `ndarray`s, where each *row* corresponds to each point). Note that in your case of only 1 input point, you don't need to transpose: your single point should be a single row vector, and `np.array([2.43‌​9, 0.097])` is equivalent to that (and actually [transposing that 1d array won't do anything](http://stackoverflow.com/a/5954747/5067311)). – Andras Deak -- Слава Україні Jan 16 '16 at 10:44