[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()
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[josef...@gmail.com]

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

[Kyle Parfrey]

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()

[josef...@gmail.com]

On Fri, Aug 20, 2010 at 2:17 PM, Kyle Parfrey <ky...@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

[Pauli Virtanen]

Fri, 20 Aug 2010 11:52:49 -0400, Kyle Parfrey wrote:
> 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):

Consider using RectBivariateSpline, if your data indeed is on a regular
grid. This way FITPACK might choose better knots...

If your data, furthermore, is given on a homogeneous grid, another
alternative is scipy.ndimage.map_coordinates. Using that requires some
thinking, though.

--
Pauli Virtanen

[Stéfan van der Walt]

On 21 August 2010 20:44, Pauli Virtanen <p...@iki.fi> wrote:
> If your data, furthermore, is given on a homogeneous grid, another
> alternative is scipy.ndimage.map_coordinates. Using that requires some
> thinking, though.

We use this routine in scikits.image.transform.homography, which
serves as a good example on how to setup the coordinates:

http://github.com/stefanv/scikits.image/blob/master/scikits/image/transform/project.py

Regards
Stéfan

[denis]

Kyle, folks,

for nonuniform 2d interpolation I'd recommend
matplotlib.mlab.griddata,
which does a Delaunay triangulation then Natural-neighbor;
see http://scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data
and http://en.wikipedia.org/wiki/Natural_neighbor .

To interpolate to arbitrary points with Delaunay,
use the little class Triinterpolate in the scipy-user thread
"Interpolation in 3D with interp2d " Aug 17.

griddata returns a masked array, masked outside the convex hull of the
input x y --
if that's a problem, you can extrapolate to the corners.

I've found fitpack BivariateSpline to be *very* sensitive to s=
(which is a sum of squares -- print sqrt( fit.get_residual() / N ),
also print once fit.get_knots() ).
What you really want to specify for BivariateSpline is the number of
knots --
but you can't, afaik.

cheers
-- denis

On Aug 20, 5:52 pm, Kyle Parfrey <kparf...@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):

[denis]

Kyle,
it looks as though there's a mixup here between Cartesion X,Y and
polar Xnew,Ynew.
griddata() will catch this to some extent with masked array Znew
but the the Fitpack routines will extrapolate *without warning*.
See the too-long notes under
http://stackoverflow.com/questions/3526514/problem-with-2d-interpolation-in-scipy-non-rectangular-grid

cheers
-- denis

[Kyle Parfrey]

Thanks everyone for all your replies, especially to Denis for those
detailed notes on stackoverflow.

I've managed to get it to work well with
interpolate.RectBivariateSpline (as suggested above) on a rectangular
(r, theta) grid, with quintic splines in theta and cubic in r, and no
smoothing (s=0). I don't know why that particular combination works so
well (quintic in both directions requires some smoothing, and the
amount needed depends on the dataset, not just the grid, which is
really bad for what I'm trying to do), I think I was being careful
enough to always convert to Cartesian coordinates (both X, Y and Xnew,
Ynew) etc, but the routines just seem to do a much better job on
rectangular grids.

If I run into problems with this later I'll try griddata() as described.

Thanks again,
Kyle