0

The following code is just to understand the context. My question doesn't require much understanding of this. It need a simple translation of one line of MATLAB code in to Python.

us = np.linspace(-(1023)/2,(1023)/2,1024)
vs = np.linspace(-(1023)/2,(1023)/2,1024)
uu,vv = np.meshgrid(-us,vs)
pu = ((((rx*SDD)/(ry+SOD))+us[0])/(-du))+1

xs = np.linspace(-(360-1)/2,(nx-1)/2,360)
ys = np.linspace(-(360-1)/2,(ny-1)/2,360)
zs = np.linspace(-(360-1)/2,(ny-1)/2,360)

xx,yy = np.meshgrid(xs,ys)

angle_rad = np.linspace(0,359,360)
angle_rad = angle_rad*np.pi/180

for i in range(0,360) :
    vol = np.zeros((360,360,360))
    rx = xx*np.cos(angle_rad[i]-np.pi/2) + yy*np.sin(angle_rad[i]-np.pi/2)
    ry = -xx*np.sin(angle_rad[i]-np.pi/2) + yy*np.cos(angle_rad[i]-np.pi/2)
    pu = ((((rx*370)/(ry+9))+us[0])/(-51.2/1024))+1

    for iz in range(0,360) :
        pv = ((((zs[iz]*370)/(ry+9))-vs[0])/(51.2/1024)) +1

So after this step the code should do interpolation and in MATLAB it's like this:

vol(:,:,iz) = interp2(uu,vv ,proj',pu,pv,'linear'); This is in MATLAB

My proj, uu and vv are (1024,1024) and pu, pv are (360,360). I need to convert the above line to Python. I tried using scipy.interpolate but it gives the following errors on trying these :

vol[:,:,iz] = Ratio*(interp2d(uu,vv,proj,pu,pv,'cubic'))

TypeError: unhashable type: 'numpy.ndarray'

vol[:,:,iz] = Ratio*(interp2d(uu,vv,proj,'cubic'))

OverflowError: Too many data points to interpolate

vol[:,:,iz] = Ratio*(interp2d(proj,pu,pv,'cubic'))

ValueError: x and y must have equal lengths for non rectangular grid

vol[:,:,iz] = Ratio*(interp2d(pu,pv,proj,'cubic'))

ValueError: Invalid length for input z for non rectangular grid

I have read all the scipy.interpolate documentations and none seemed to help. Could anyone figure out what's wrong?

Prakriti Kumari
  • 49
  • 1
  • 12
  • 1
    Read http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html and look at the example at the bottom of the page. Rather than just reading it, open a Python console and repeat each of the steps and try printing the arrays. – pv. Mar 12 '16 at 13:06

1 Answers1

1

The problem on a wide scale is that you're looking at the documentation but you're not trying to understand it. If you're using a specific function interp2d in the module scipy.interpolate, then look at the function's documentation, as @pv also suggested in a comment. The fact that you're throwing the arguments all around the place clearly demonstrates that you're trying to use a function based on guesswork. It won't work: a function was implemented with a given syntax, and it will only work like that.

So look at the function's signature:

class scipy.interpolate.interp2d(x, y, z, kind='linear', copy=True, bounds_error=False, fill_value=nan)

The meaning of the parameters is explained afterwards. You can clearly see that there are 3 mandatory parameters: x, y, z. No other array-valued inputs are allowed This is because interp2d only constructs an interpolating function, which you should then use to compute the interpolated values on a mesh (unlike MATLAB, where interp2d gives you the interpolated values directly). So you can call

myfun = interp2(uu,vv,proj,'linear')

to get an interpolating function. You can then substitute at the given values, but note that the function myfun will expect 1d input, and it will construct a mesh internally. So assuming that your output mesh is constructed as

puu,puv = np.meshgrid(puu_vec,puv_vec)

(which is probably not the case, but I'll get back to this later), you need

vol[:,:,iz] = Ratio*myfun(puu_vec,puv_vec)

to obtain the output you need.

However, there are some important points you should note.

  1. In MATLAB you have interp2d(uu,vv,proj',...), but make sure that for scipy the elements of uu, vv, and proj are in the same order. To be on the safe side: in case of an asymmetric mesh size, the shape of uu, vv and proj should all be the same.
  2. In MATLAB you're using 'linear' interpolation, while in python you're using 'cubic'. I'm not sure this is really what you want, if you're porting your code to a new language.
  3. It seems to me that your output mesh is not defined by a rectangular grid as if from meshgrid, which suggests that interp2d might not be suitable for your case. Anyway, I've had some odd experiences with interp2d, and I wouldn't trust it. So if it's not suitable for your expected output, I'd strongly suggest using scipy.interpolate.griddata instead. This function gives you interpolated points directly: I suggest that you try to figure out its use based on the manual:) My linked answer can also help. You can set the kind of interpolation in the same way, but your output can be any set of scattered points if you like.
Community
  • 1
  • 1
  • I tried this. f = interp2d(uu,vv,proj,'cubic') vol[:,:,iz] = Ratio*f(pu,pv) But it says OverflowError: Too many data points to interpolate. My pu and pv are (360,360) and if the inputs are supposed to be 1-D I tried pu=pu.reshape(pu.size) and same for pv to make it 1-D. But it still gives the same overflow error. – Prakriti Kumari Mar 13 '16 at 08:42
  • @PrakritiKumari read my example code carefully. It's not the *shape* that is the problem, but the structure of the points themselves. If your input points would be `[0,0]`, `[0,1]`, `[1,0]`, `[1,1]`, then your `pu` and `pv` are something like `[[0,0],[1,1]]` and `[[0,1],[0,1]]` (or the other way around). If you want to interpolate over these 2x2 points, you have to give only the two 2-element vectors to the interpolator: `f([0,1],[0,1])`. So the same vectors go into `f()` that can be used to generate your mesh with `meshgrid`. – Andras Deak -- Слава Україні Mar 13 '16 at 10:46
  • This is why I'm saying that if your points are not on a grid, you can't use `interp2d`. Or at least you're way better off with `griddata`. Do you understand now? – Andras Deak -- Слава Україні Mar 13 '16 at 10:47