user2063407
user2063407

Reputation: 41

How can I interpolate array from spherical to cartesian coordinates with Python?

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 density2 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.

Upvotes: 4

Views: 3686

Answers (1)

Dave
Dave

Reputation: 8090

Something like this should work:

import scipy.interpolate
rflat=scipy.array( r_coord.flat )
tflat=scipy.array( theta_coord.flat )
pflat=scipy.array( phi_coord.flat )
coordpoints=scipy.concatenate( [ rflat[:, scipy.newaxis], tflat[:,scipy.newaxis], pflat[:,scipy.newaxis] ] , axis=1 )
rtpinterpolator=scipy.interpolate.linearNDInterpolate( coordppoints, density.flat )    

def xyz2rtp( x,y,z):
     r=scipy.sqrt( x**2+y**2+z**2)
     t=scipy.acos( z/r )
     p=scipy.atan2( y, x )
     return (r,t,p)

# now you can get the interpolated value for any (x,y,z) coordinate you want.
val=rtpinterpolator( xyz2rtp( x,y,z) )

Key points:

  • Use the existing scipy multi-dimensional interpolation,
  • convert the xyz coordinates to rtp when you pass it in to the interpolator.

Upvotes: 2

Related Questions