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

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Joseph Smidt-3
I have an array of density values in spherical coordinates. More specifically I have an array called density with shape (180,200,200). I also have an array called r_coord, theta_coord and phi_coord also with shape (180,200,200) being the spherical coordinates for the density array.

I would like to map this density to cartesian coordinates using python. I will need therefore a new array density_prime which is interpolated over cartesian coordinates x_coord, y_coord and z_coord. I found scipy.ndimage.interpolation.map_coordinates which looks promising but I can't figure out how to get it to work.

Any help would be appreciated. Thanks.

--
------------------------------------------------------------------------
Joseph Smidt <[hidden email]>

Theoretical Division
P.O. Box 1663, Mail Stop B283
Los Alamos, NM 87545
Office: 505-665-9752
Fax:    505-667-1931

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: How can I interpolate array from spherical to cartesian coordinates?

Charles R Harris


On Mon, Feb 11, 2013 at 8:48 PM, Joseph Smidt <[hidden email]> wrote:
I have an array of density values in spherical coordinates. More specifically I have an array called density with shape (180,200,200). I also have an array called r_coord, theta_coord and phi_coord also with shape (180,200,200) being the spherical coordinates for the density array.

I would like to map this density to cartesian coordinates using python. I will need therefore a new array density_prime which is interpolated over cartesian coordinates x_coord, y_coord and z_coord. I found scipy.ndimage.interpolation.map_coordinates which looks promising but I can't figure out how to get it to work.


I'm not clear on what you are trying to do, but I'm guessing you have sample points on a sphere and you want to find interpolated values at other points on the sphere, the cartesian coordinates being a means rather than an end. Is that the case? If not, can you be more explicit.

Chuck

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: How can I interpolate array from spherical to cartesian coordinates?

Joseph Smidt-3
Chuck and everyone,

   Okay, I will give a more specific example.  Consider the following script:

from pylab import *

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

r = r/20.0                     # r goes from 0 to 10.
theta = theta/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)

# Plot density.  Doesn't look spherical because 
# not in cartesian coordinates.
imshow(squeeze(f[:,:,1]))


   Okay, so density is defined in terms of spherical coordinates.  I would like a function that transforms density to density_prime to cartesian coordinate arrays x, y, and z such that the r = 0 line gets mapped to x = 0, y = 0 z = 0.  The r = 1 line gets mapped to x = sin(theta)*cos(phi), y = sin(theta)*sin(phi), z = cos(phi), the r = 2 line gets mapped to... etc.  So that when I plot density_prime it looks like a nice spherical function peaking at x = 0, y = 0 z = 0 and getting small for x, y, and z large.

  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.





On Mon, Feb 11, 2013 at 11:07 PM, Charles R Harris <[hidden email]> wrote:


On Mon, Feb 11, 2013 at 8:48 PM, Joseph Smidt <[hidden email]> wrote:
I have an array of density values in spherical coordinates. More specifically I have an array called density with shape (180,200,200). I also have an array called r_coord, theta_coord and phi_coord also with shape (180,200,200) being the spherical coordinates for the density array.

I would like to map this density to cartesian coordinates using python. I will need therefore a new array density_prime which is interpolated over cartesian coordinates x_coord, y_coord and z_coord. I found scipy.ndimage.interpolation.map_coordinates which looks promising but I can't figure out how to get it to work.


I'm not clear on what you are trying to do, but I'm guessing you have sample points on a sphere and you want to find interpolated values at other points on the sphere, the cartesian coordinates being a means rather than an end. Is that the case? If not, can you be more explicit.

Chuck

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user




--
------------------------------------------------------------------------
Joseph Smidt <[hidden email]>

Theoretical Division
P.O. Box 1663, Mail Stop B283
Los Alamos, NM 87545
Office: 505-665-9752
Fax:    505-667-1931

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: How can I interpolate array from spherical to cartesian coordinates?

Jerome Kieffer
On Mon, 11 Feb 2013 23:30:11 -0700
Joseph Smidt <[hidden email]> 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
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: How can I interpolate array from spherical to cartesian coordinates?

Joseph Smidt-3
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 <[hidden email]> wrote:
On Mon, 11 Feb 2013 23:30:11 -0700
Joseph Smidt <[hidden email]> 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 <a href="tel:%2B33%20476%20882%20445" value="+33476882445">+33 476 882 445



--
------------------------------------------------------------------------
Joseph Smidt <[hidden email]>

Theoretical Division
P.O. Box 1663, Mail Stop B283
Los Alamos, NM 87545
Office: 505-665-9752
Fax:    505-667-1931

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user