[SciPy-User] Problems with 2D interpolation of data on polar grid

[Kyle Parfrey]

Hi all,

I've been having some trouble doing 2D interpolation with both
interp2d and bisplrep/bisplev. My data is on a spherical polar (r,
theta) grid, and I'm trying to interpolate functions similar to the
vector components of a dipole field. The output looks very wrong, and
I keep getting warnings like (from interp2d):

Warning:     No more knots can be added because the number of B-spline
coefficients
    already exceeds the number of data points m. Probably causes: either
    s or m too small. (fp>s)
	kx,ky=1,1 nx,ny=16,9 m=90 fp=0.000000 s=0.000000

or, from bisplrep:

/usr/lib/python2.6/dist-packages/scipy/interpolate/fitpack.py:763:
DeprecationWarning: integer argument expected, got float
  tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)

One thing: I'm not trying to regrid---I'll need to be able to
interpolate to any arbitrary point, since I'll be using this to
calculate field lines. I've pasted some code below that should give
some idea of the problem. I hope I haven't done something obviously
stupid!

Thanks in advance,
Kyle

----------------------------------------------------

import numpy as np
from scipy import interpolate
from math import *
import matplotlib.pyplot as plt

### Make polar grid ###
rvec = np.arange(1.0, 11.0, 1.0)
tvec = np.arange(pi/10.0, pi, pi/10.0)
Nr = len(rvec)
Nt = len(tvec)
X = np.empty([Nr,Nt])
Y = np.empty([Nr,Nt])
Z = np.empty([Nr,Nt])

for i in range(Nr):
  for j in range(Nt):
    r = rvec[i]
    t = tvec[j]
    X[i,j] = r*sin(t)
    Y[i,j] = r*cos(t)
    Z[i,j] = cos(t)/pow(r,3) # cos(theta)/r^3: Br of dipole


### Do the interpolation ###
#interp_poly = interpolate.interp2d(X,Y,Z, kind='linear')
tck = interpolate.bisplrep(X,Y,Z, kx=3, ky=3)


### interpolate onto new grid ###
rnew = np.arange(1.0, 10.1, 0.1)
tnew = np.arange(pi/100.0, pi, pi/100.0)
Nr2 = len(rnew)
Nt2 = len(tnew)
X2 = np.empty([Nr2,Nt2])
Y2 = np.empty([Nr2,Nt2])
Z2 = np.empty([Nr2,Nt2])
for i in range(Nr2):
  for j in range(Nt2):
    r = rnew[i]
    t = tnew[j]
    X2[i,j] = r*sin(t)
    Y2[i,j] = r*cos(t)
    #Z2[i,j] = interp_poly(X2[i,j], Y2[i,j])
    Z2[i,j] = interpolate.bisplev(X2[i,j], Y2[i,j], tck)


### Pseudocolour plot ###
fig = plt.figure()
fig.add_subplot(111, aspect='equal')
plt.pcolor(X2,Y2,Z2)
plt.colorbar()
plt.show()