[SciPy-User] How can I interpolate array from spherical to cartesian coordinates?

[Joseph Smidt]

Hey everyone,

      I found a solution to this problem so I decided to post it here for
posterity's sake.  The code below seems to do the trick:

from pylab import *
from scipy.interpolate import interp1d
from scipy.ndimage import map_coordinates
import scitools

def spherical2cartesian(r, th, phi, grid, x, y, z, order=3):

    # Build relationship between Cartesian and spherical coordinates.
    X, Y, Z = scitools.numpytools.meshgrid(x,y,z)
    new_r = np.sqrt(X*X+Y*Y+Z*Z)
    new_th = np.arccos(Z/new_r)
    new_phi = np.arctan2(Y, X)

    # Find these values for the input grid
    ir = interp1d(r, np.arange(len(r)), bounds_error=False)
    ith = interp1d(th, np.arange(len(th)))
    iphi = interp1d(phi, np.arange(len(phi)))
    new_ir = ir(new_r.ravel())

    new_ith = ith(new_th.ravel())
    new_iphi = iphi(new_phi.ravel())

    new_ir[new_r.ravel() > r.max()] = len(r)-1
    new_ir[new_r.ravel() < r.min()] = 0

    # Interpolate to Cartesian coordinates.
    return map_coordinates(grid, np.array([new_ir, new_ith, new_iphi]),
                            order=order).reshape(new_r.shape)

# Build 3D arrays for spherical coordinates.
r, th, phi = mgrid[0:201,0:201,0:201]

r = r/20.0               # r goes from 0 to 10.
th = th/200.0*pi         # Theta goes from 0 to pi
phi = phi/200.0*2*pi     # Phi goes from 0 to 2pi

# Density is spherically symmetric.  Only depends on r.
density = exp(-r**2/20.0)

# Build ranges for function
r = linspace(0,200,200)
th = linspace(0,np.pi,200)
phi = linspace(-np.pi,np.pi,200)
x = linspace(-200,200,200)
y = linspace(-200,200,200)
z = linspace(-200,200,200)

# Map to Cartesian coordinates.
density = spherical2cartesian(r, th, phi, density, x, y, z, order=3)

# Plot density, now in spherical coordinates.
# not in cartesian coordinates.
figure()
imshow(squeeze(density[:,:,100]))
show()



On Tue, Feb 12, 2013 at 1:02 AM, Jerome Kieffer <Jerome.Kieffer at esrf.fr>wrote:

> On Mon, 11 Feb 2013 23:30:11 -0700
> Joseph Smidt <josephsmidt at gmail.com> wrote:
>
>
> >   If anyone knows how I could do such a transformation to get
> density_prime
> > with scipy.ndimage.interpolation.map_coordinates or any other
> interpolator
> > for N-dim data I would appreciate it.
>
> I never did it in 3D but you need the inverse transformation for
> map_coordinate:
> r, theta, phi -> x, y , z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi),
> r*cos(phi)
> I think that's all.
>
> Cheers,
> --
> Jérôme Kieffer
> On-Line Data analysis / Software Group
> ISDD / ESRF
> tel +33 476 882 445
>



-- 
------------------------------------------------------------------------
Joseph Smidt <josephsmidt at gmail.com>

Theoretical Division
P.O. Box 1663, Mail Stop B283
Los Alamos, NM 87545
Office: 505-665-9752
Fax:    505-667-1931
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130212/6f9e8f12/attachment.html>