Faster numpy cartesian to spherical coordinate conversion?

Question

I have an array of 3 million data points from a 3-axiz accellerometer (XYZ), and I want to add 3 columns to the array containing the equivalent spherical coordinates (r, theta, phi). The following code works, but seems way too slow. How can I do better?

import numpy as np import math as m  def cart2sph(x,y,z):     XsqPlusYsq = x**2 + y**2     r = m.sqrt(XsqPlusYsq + z**2)               # r     elev = m.atan2(z,m.sqrt(XsqPlusYsq))     # theta     az = m.atan2(y,x)                           # phi     return r, elev, az  def cart2sphA(pts):     return np.array([cart2sph(x,y,z) for x,y,z in pts])  def appendSpherical(xyz):     np.hstack((xyz, cart2sphA(xyz))) 

Answer 1

This is similar to Justin Peel's answer, but using just numpy and taking advantage of its built-in vectorization:

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 

Note that, as suggested in the comments, I've changed the definition of elevation angle from your original function. On my machine, testing with pts = np.random.rand(3000000, 3), the time went from 76 seconds to 3.3 seconds. I don't have Cython so I wasn't able to compare the timing with that solution.