如何用Python将数组从球坐标插值到笛卡尔坐标？

Question

我有一个表示密度值的数组，这些值是用球坐标表示的。更具体地说，我有一个叫做 density 的数组，它的形状是 (180,200,200)。我还有三个数组，分别叫做 r_coord、theta_coord 和 phi_coord，它们的形状也是 (180,200,200)，用来表示这个密度数组的球坐标。

我想用 Python 把这个密度映射到笛卡尔坐标系中。因此，我需要一个新的数组 density2，它是在笛卡尔坐标 x_coord、y_coord 和 z_coord 上进行插值后的结果。我发现了 scipy.ndimage.interpolation.map_coordinates 这个函数，看起来很有用，但我不知道怎么用它。

如果有人能帮忙就太好了。谢谢！

Answer 1

像这样应该可以工作：

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

关键点：

• 使用现有的scipy多维插值功能，
• 在传入插值器时，把xyz坐标转换成rtp格式。