Most efficient way to generate a large array of (x,y,z) coordinates

user_1663725 0
Viewed 78

I'm generating coordinates of a bi-cone model in spherical coordinates. I'm using a series of nested for loops, as show here:

theta_in  = 30.0 * np.pi/180.0
theta_out = 60.0 * np.pi/180.0
phi       = 2*np.pi # rotation
R = 1.0
sampling = 100

theta = np.linspace(theta_in,theta_out,sampling)
phi   = np.linspace(0,phi,sampling)
r     = np.linspace(-R,R,sampling)

x = []
y = []
z = []
for ri in r:
    for pi in phi:
         for ti in theta:
            xi = ri*np.cos(pi)*np.sin(ti)
            yi = ri*np.sin(pi)*np.sin(ti)
            zi = ri*np.cos(ti)
            x.append(xi)
            y.append(yi)
            z.append(zi)

This generates the following desired output:

enter image description here

The use of these nested for loops to generate coordinate lists is not very efficient. My question is: what would be the most efficient way to do this while possibly avoiding the use of nested for loops?

numpy arrays python
2 Answers

Create open grids off those inputs and then perform the same operations -

RI,PI,TI = np.ix_(r,phi,theta) # get open grids           
X = RI*np.cos(PI)*np.sin(TI)
Y = RI*np.sin(PI)*np.sin(TI)
Z = np.repeat(RI*np.cos(TI),sampling,axis=1)

Alternative #1 : Those open grids could also be constructed with explicit axis addition, like so -

RI,PI,TI = r[:,None,None], phi[:,None], theta

Alternative #2 : We can pre-compute np.sin(TI) and re-use it at two steps.

Alternative #3 : Also, for a read-only broadcatsed version of Z, we can use np.broadcast_to -

Z = np.broadcast_to(RI*np.cos(TI),(sampling,sampling,sampling))

Alternative #4 : Leverage multi-cores with numexpr module -

import numexpr as ne

X = ne.evaluate('RI*cos(PI)*sin(TI)')
# similarly Y and Z

Note that the outputs would be 3D arrays. So, to get equivalent ones, you might want to flatten those. For the same, we can use .ravel() on the outputs.

user_32938810 
theta_in  = 30.0 * np.pi/180.0
theta_out = 60.0 * np.pi/180.0
phi       = 2*np.pi # rotation
R = 1.0
sampling = 100

theta = np.linspace(theta_in,theta_out,sampling)
phi   = np.linspace(0,phi,sampling)
r     = np.linspace(-R,R,sampling)

x=list(map(lambda ti,pi,ri:ri*np.cos(pi)*np.sin(ti),theta,phi,r))
y=list(map(lambda ti,pi,ri:ri*np.sin(pi)*np.sin(ti),theta,phi,r))
z=list(map(lambda ti,ri:ri*np.cos(ti),theta,r))
user_121600790
Related