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

Kyle Parfrey kyle@astro.columbia....
Fri Aug 20 13:17:29 CDT 2010


The interpolated field doesn't look as crazy in the centre of the
grid, but I think it's still wrong. Have a go with the code below, in
which I've fixed a few little things, and
now you see the original field on the left, interpolated on the right.
(I'm having trouble forcing the colour ranges to be the same in both
plots.)

Is the "hole" in the middle of the grid causing problems?

Thanks,
Kyle


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

### Make polar grid ###
rvec = np.linspace(1.0, 5.0, num=50)
tvec = np.linspace(0, pi, num=50)
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.linspace(1.0, 5.0, num=100)
tnew = np.linspace(0, pi, num=100)
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)


### Plot pseudocolour plot ###
fig = plt.figure()
fig.add_subplot(121, aspect='equal')
plt.pcolor(X,Y,Z)
plt.colorbar()

fig.add_subplot(122, aspect='equal')
plt.pcolor(X2,Y2,Z2)
plt.colorbar()
plt.show()






On 20 August 2010 12:58,  <josef.pktd@gmail.com> wrote:
> On Fri, Aug 20, 2010 at 11:52 AM, Kyle Parfrey <kparfrey@gmail.com> wrote:
>> 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()
>
>
> I'm not a graphical person, but if I reduce the radius
>
> ### interpolate onto new grid ###
> rnew = np.arange(1.0, 10.1, 0.1)[:20]
>
> the plot seems to look reasonable.
>
> plotting the original X,Y,Z also has color variation only close to the
> origin. contour plot also shows "action" only around the origin.
> looks like most of the area is flat close to zero.
>
> I don't see a problem with the spline interpolation itself.
>
> Josef
>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list