Более быстрое преобразование декартовых координат в сферические координаты?

Это похоже на ответ Justin Peel, но используя только numpy и используя встроенную векторизацию:

import numpy as np

def appendSpherical_np(xyz):
    ptsnew = np.hstack((xyz, np.zeros(xyz.shape)))
    xy = xyz[:,0]**2 + xyz[:,1]**2
    ptsnew[:,3] = np.sqrt(xy + xyz[:,2]**2)
    ptsnew[:,4] = np.arctan2(np.sqrt(xy), xyz[:,2]) # for elevation angle defined from Z-axis down
    #ptsnew[:,4] = np.arctan2(xyz[:,2], np.sqrt(xy)) # for elevation angle defined from XY-plane up
    ptsnew[:,5] = np.arctan2(xyz[:,1], xyz[:,0])
    return ptsnew

Обратите внимание, что, как было предложено в комментариях, я изменил определение угла возвышения по сравнению с исходной функцией. На моей машине при тестировании с pts = np.random.rand(3000000, 3) время увеличилось с 76 до 3,3 секунды. У меня нет Cython, поэтому я не смог сравнить время с этим решением.

Вот быстрый код Cython, который я написал для этого:

cdef extern from "math.h":
    long double sqrt(long double xx)
    long double atan2(long double a, double b)

import numpy as np
cimport numpy as np
cimport cython

ctypedef np.float64_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
def appendSpherical(np.ndarray[DTYPE_t,ndim=2] xyz):
    cdef np.ndarray[DTYPE_t,ndim=2] pts = np.empty((xyz.shape[0],6))
    cdef long double XsqPlusYsq
    for i in xrange(xyz.shape[0]):
        pts[i,0] = xyz[i,0]
        pts[i,1] = xyz[i,1]
        pts[i,2] = xyz[i,2]
        XsqPlusYsq = xyz[i,0]**2 + xyz[i,1]**2
        pts[i,3] = sqrt(XsqPlusYsq + xyz[i,2]**2)
        pts[i,4] = atan2(xyz[i,2],sqrt(XsqPlusYsq))
        pts[i,5] = atan2(xyz[i,1],xyz[i,0])
    return pts

Мне удалось сократить время с 62,4 секунды до 1,22 секунды, используя 3 000 000 баллов. Это не так уж и плохо. Я уверен, что есть и другие улучшения, которые можно сделать.

! Во всем вышеприведенном коде по-прежнему есть ошибка... и это лучший результат Google.. TLDR: я тестировал это с VPython, использование atan2 для тета (уровень) неверно, используйте acos! Правильно для фи (азим). Я рекомендую функцию acos sympy1.0 (она даже не жалуется на acos(z/r) с r = 0).

http://mathworld.wolfram.com/SphericalCoordinates.html

Если мы преобразуем это в физическую систему (r, тета, фи) = (r, высота над уровнем моря, азимут), мы получим:

r = sqrt(x*x + y*y + z*z)
phi = atan2(y,x)
theta = acos(z,r)

Неоптимизированный, но правильный код для правой системы физики:

from sympy import *
def asCartesian(rthetaphi):
    #takes list rthetaphi (single coord)
    r       = rthetaphi[0]
    theta   = rthetaphi[1]* pi/180 # to radian
    phi     = rthetaphi[2]* pi/180
    x = r * sin( theta ) * cos( phi )
    y = r * sin( theta ) * sin( phi )
    z = r * cos( theta )
    return [x,y,z]

def asSpherical(xyz):
    #takes list xyz (single coord)
    x       = xyz[0]
    y       = xyz[1]
    z       = xyz[2]
    r       =  sqrt(x*x + y*y + z*z)
    theta   =  acos(z/r)*180/ pi #to degrees
    phi     =  atan2(y,x)*180/ pi
    return [r,theta,phi]

вы можете проверить это самостоятельно с помощью такой функции, как:

test = asCartesian(asSpherical([-2.13091326,-0.0058279,0.83697319]))

некоторые другие тестовые данные для некоторых квадрантов:

[[ 0.          0.          0.        ]
 [-2.13091326 -0.0058279   0.83697319]
 [ 1.82172775  1.15959835  1.09232283]
 [ 1.47554111 -0.14483833 -1.80804324]
 [-1.13940573 -1.45129967 -1.30132008]
 [ 0.33530045 -1.47780466  1.6384716 ]
 [-0.51094007  1.80408573 -2.12652707]]

Я использовал VPython дополнительно, чтобы легко визуализировать векторы:

test   = v.arrow(pos = (0,0,0), axis = vis_ori_ALA , shaftwidth=0.05, color=v.color.red)

Чтобы завершить предыдущие ответы, вот реализация Numexpr (с возможным откатом к Numpy),

import numpy as np
from numpy import arctan2, sqrt
import numexpr as ne

def cart2sph(x,y,z, ceval=ne.evaluate):
    """ x, y, z :  ndarray coordinates
        ceval: backend to use: 
              - eval :  pure Numpy
              - numexpr.evaluate:  Numexpr """
    azimuth = ceval('arctan2(y,x)')
    xy2 = ceval('x**2 + y**2')
    elevation = ceval('arctan2(z, sqrt(xy2))')
    r = eval('sqrt(xy2 + z**2)')
    return azimuth, elevation, r

Для больших размеров массива это позволяет увеличить скорость в 2 раза по сравнению с чистой реализацией Numpy и будет сопоставимо со скоростью C или Cython. Текущее решение numpy (при использовании с аргументом ceval=eval) также на 25% быстрее, чем функция appendSpherical_np в ответе @mtrw для больших размеров массива,

In [1]: xyz = np.random.rand(3000000,3)
   ...: x,y,z = xyz.T
In [2]: %timeit -n 1 appendSpherical_np(xyz)
1 loops, best of 3: 397 ms per loop
In [3]: %timeit -n 1 cart2sph(x,y,z, ceval=eval)
1 loops, best of 3: 280 ms per loop
In [4]: %timeit -n 1 cart2sph(x,y,z, ceval=ne.evaluate)
1 loops, best of 3: 145 ms per loop

хотя для меньших размеров appendSpherical_np на самом деле быстрее,

In [5]: xyz = np.random.rand(3000,3)
...: x,y,z = xyz.T
In [6]: %timeit -n 1 appendSpherical_np(xyz)
1 loops, best of 3: 206 µs per loop
In [7]: %timeit -n 1 cart2sph(x,y,z, ceval=eval)
1 loops, best of 3: 261 µs per loop
In [8]: %timeit -n 1 cart2sph(x,y,z, ceval=ne.evaluate)
1 loops, best of 3: 271 µs per loop

Octave имеет некоторые встроенные функции для преобразования координат, к которым можно получить доступ с помощью пакета oct2py для преобразования массивов numpy в декартовых координатах в сферические или полярные координаты (и обратно):

from oct2py import octave
xyz = np.random.rand(3000000,3)
%timeit thetaphir = octave.cart2sph(xyz)

724 ms ± 206 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)