Problems interpolating and evaluating numpy array at arbitrary points with Scipy

Question

I am trying to replicate some of the functionality of Matlab's interp2. I know somewhat similar questions have been asked before, but none apply to my specific case.

I have a distance map (available at this Google drive location): https://drive.google.com/open?id=0B6acq_amk5e3X0Q5UG1ya1VhSlE&authuser=0

Values are normalized in the range 0-1. Size is 200 rows by 300 columns.

I can load it up with this code snippet:

import numpy as np
dstnc1=np.load('dstnc.npy')

Coordinates are defined by the next snippet:

xmin = 0.
xmax = 9000.
ymin = 0.
ymax = 6000.
r1,c1 = dstnc1.shape
x = np.linspace(xmin,xmax,c1)
y = np.linspace(ymin, ymax,r1)

I have three map points defined by vectors xnew1, ynew1 with this snippet:

xnew1=[3700.540199,3845.940199,3983.240199]
ynew1=[1782.8611,1769.862,1694.862]

I check their location with respect to the distance map with this:

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(20, 16))
ax = fig.add_subplot(1, 1, 1)
plt.imshow(dstnc1, cmap=my_cmap_r,vmin=0,vmax=0.3,
           extent=[0, 9000, 0, 6000], origin='upper')
plt.scatter(xnew1, ynew1, s=50, linewidths=0.15)
plt.show()

They plot in the correct location. Now I would like to extract the distance value at those three points. I tried first interp2d.

from scipy.interpolate import interp2d
x1 = np.linspace(xmin,xmax,c1)
y1 = np.linspace(ymin,ymax,r1)
f = interp2d(x1, y1, dstnc1, kind='cubic')

but when I try to evaluate with:

test=f(xnew1,ynew1)

I get this error:

--------------------
ValueError                                Traceback (most recent call last)
<ipython-input-299-d0f42e609b23> in <module>()
----> 1 test=f(xnew1,ynew1)

C:\...\AppData\Local\Continuum\Anaconda\lib\site-packages\scipy\interpolate\interpolate.pyc
in __call__(self, x, y, dx, dy)
    270                                 (self.y_min, self.y_max)))
    271
--> 272         z = fitpack.bisplev(x, y, self.tck, dx, dy)
    273         z = atleast_2d(z)
    274         z = transpose(z)

C:\...\AppData\Local\Continuum\Anaconda\lib\site-packages\scipy\interpolate\fitpack.pyc
in bisplev(x, y, tck, dx, dy)
   1027     z,ier = _fitpack._bispev(tx,ty,c,kx,ky,x,y,dx,dy)
   1028     if ier == 10:
-> 1029         raise ValueError("Invalid input data")
   1030     if ier:
   1031         raise TypeError("An error occurred")

ValueError: Invalid input data

If I try RectBivariateSpline:

from scipy.interpolate import RectBivariateSpline
x2 = np.linspace(xmin,xmax,r1)
y2 = np.linspace(ymin,ymax,c1)
f = RectBivariateSpline(x2, y2, dstnc1)

I get this error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-302-d0f42e609b23> in <module>()
----> 1 test=f(xnew1,ynew1)

C:\...\AppData\Local\Continuum\Anaconda\lib\site-packages\scipy\interpolate\fitpack2.pyc
in __call__(self, x, y, mth, dx, dy, grid)
    643                 z,ier = dfitpack.bispev(tx,ty,c,kx,ky,x,y)
    644                 if not ier == 0:
--> 645                     raise ValueError("Error code returned by
bispev: %s" % ier)
    646         else:
    647             # standard Numpy broadcasting

ValueError: Error code returned by bispev: 10

Any suggestion as to whether I am using the wrong functions or the right function with wrong syntax, and how I may fix it is appreciated. Thank you.

UPDATE

I am running Python 2.7.9 and Scipy 0.14.0 (on Continuum Anaconda) As posted on the Scipy mailing list here the documentation seems confusing, being a mix of Scipy 0.14.0, and the next version. Can anybody suggest a workaround or the correct syntax for version 0.14.0.

UPDATE 2

tried

xnew1=np.array([3700.540199,3845.940199,3983.240199])
ynew1=np.array([1782.8611,1769.862,1694.862])

as suggested inj a comment but the error remains.

Tried your code and it worked. What scipy version do you use ? — Moritz, May 14 '15 at 23:41
It is interesting to see how colors can introduce artefacts. Try these colormaps and see the difference: `cm = plt.cm.BrBG #cm = plt.cm.gray #cm = plt.cm.cubehelix` Can you comment on that ? — Moritz, May 14 '15 at 23:45
+1 about colour ,Moritz I could not agree more. I wrote quite a bit about it http://nbviewer.ipython.org/github/mycarta/tutorials/blob/master/1408_Evaluate_and_compare_colormaps/How_to_evaluate_and_compare_colormaps.ipynb and often use one very similar to cubehelix. I noticed that with distance maps artifacts are even more pronounced, probably because it is strictly increasing data like the pyramid in my tutorial. — MyCarta, May 15 '15 at 02:17
About the xnew, ynew instead of xnew1, ynew1, I must have copied that part from the older notebook, but I am nearly certain I had the latter in the newest. I'll have to check tomorrow in the office and also check the scipy version on that PC. I will add as an edit to the question accordingly. — MyCarta, May 15 '15 at 02:25
Moritz, I tried again with xnew1, ynew1, still get the errors. I run I am running Python 2.7.9 and Scipy 0.14.0, see update. Can you think of a workaround? — MyCarta, May 15 '15 at 14:02
I use scipy 15.01 within the conda environment (continuum analytics). What if you provide xnew1 and ynew1 as numpy array and not as list ? — Moritz, May 15 '15 at 19:58
Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/77936/discussion-between-mycarta-and-moritz). — MyCarta, May 15 '15 at 22:34

score 1 · Accepted Answer · answered May 20 '15 at 21:56

This syntax worked with RectBivariateSpline

x2 = np.linspace(xmin,xmax,c1)
y2 = np.linspace(ymin,ymax,r1)

f2 = sp.interpolate.RectBivariateSpline(x2, y2, dstnc1.T,kx=1, ky=1)

I can then evaluate at new points with this:

out2 = f2.ev(xnew1,ynew1)

For interp2d I am stuck as I am not able to bypass firewall at my office to update Anaconda (Windows). I may be able at home on a Mac installation, in which case, if I get the syntax right, I will add to thsi answer.

Problems interpolating and evaluating numpy array at arbitrary points with Scipy

UPDATE

UPDATE 2

1 Answers1