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

[josef.pktd at gmail.com]

On Fri, Aug 20, 2010 at 2:17 PM, Kyle Parfrey <kyle at astro.columbia.edu> wrote:
> 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?

maybe the radial placement of the points ? I don't know.

adding a small positive s seems to help in the graph  (I tried s=
0.001 to s=0.1)

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

looks much closer to the original but I didn't check what error you
get at the original points.

Josef


>
> 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 at gmail.com> wrote:
>> On Fri, Aug 20, 2010 at 11:52 AM, Kyle Parfrey <kparfrey at 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 at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>