2

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.

10
  • 2
    Tried your code and it worked. What scipy version do you use ? Commented May 14, 2015 at 23:41
  • 2
    You did use xǹew . Did you tried xnew1 ? Commented May 14, 2015 at 23:43
  • 1
    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 ? Commented May 14, 2015 at 23:45
  • +1 about colour ,Moritz I could not agree more. I wrote quite a bit about it nbviewer.ipython.org/github/mycarta/tutorials/blob/master/… 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. Commented May 15, 2015 at 2:17
  • 1
    what if you do conda update scipy ? I think it is a bug Commented May 15, 2015 at 21:51

1 Answer 1

1

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.

Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.