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 |
On Mon, Feb 11, 2013 at 8:48 PM, Joseph Smidt <[hidden email]> wrote:
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 |
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:
------------------------------------------------------------------------ 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 |
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 |
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:
------------------------------------------------------------------------ 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 |
Free forum by Nabble | Edit this page |