Not sure if this question has been asked before–I looked through similar examples and they weren't exactly what I need to do.

I have an array of positions (shape = (8855470, 3)) in a cube with physical coordinates in between 0 and 787.5. These positions represent point masses in some space. Here's a look at the first three entries of this array:

array([[224.90635586, 720.494766  ,  19.40263367],
   [491.25279546,  41.26026654,   7.35436416],
   [407.70436788, 340.32618713, 328.88192913]])

I want to split this giant cube into a number of smaller cubes. For example, if I wanted to split it on each side into 10 cubes, making 1,000 subcubes total, then each subcube would contain only the points that have positions within that subcube. I have been experimenting with np.meshgrid to create the 3D grid necessary to conditionally apportion the appropriate entries of the positions array to subcubes:

split = np.arange(0.,(787.5 787.5/10.),step=787.5/10.)
xg,yg,zg = np.meshgrid(split,split,split,indexing='ij')

But I'm not sure if this is the way to go about this. Let me know if this question is too vague or if you need any additional information.

CodePudding user response：

For sake of problem I will work with toy data. I think you're near with the meshgrid. Here's a propossal

1. Create grid but with points until 757.5 not included, with values as you did in arange.
2. Reshape then to have a 1d_array. for in arrays zip to get masks with the cube shape.
3. create a list to save all subcube points.

    import numpy as np
    data = np.random.randint(0,787,( 10000,3))
    
    start = 0
    end = 787.5
    step = (end-start)/10
    split = np.arange(start,end,step)
    
    xg,yg,zg = np.meshgrid(split,split,split,indexing='ij')
    
    xg = xg.reshape(-1)
    yg = yg.reshape(-1)
    zg = zg.reshape(-1)

    subcube_data = []
    for x,y,z in zip(xg,yg,zg):
        mask_x = (x<= data[:,0] ) * ( data[:,0] < x step) #data_x between start and end for this subcube
        mask_y = (y<= data[:,1] ) * ( data[:,1] < y step) #data_y between start and end for this subcube
        mask_z = (z<= data[:,2] ) * ( data[:,2] < z step) #data_z between start and end for this subcube
        mask = mask_x * mask_y * mask_z
        subcube_data.append(data[mask])

Now you will have a list with 1000 elements where each element is a sub_cube containing an Nx3 point list. If you want to recover the 3d index corresponding to every sub_cube[i] you just could do [xg[i],yg[i],zg[i]].

Last you can plot to see some of the sub_cubes and the rest of data

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#plot data as 3d scatter border black
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

#plot subcubes 0 1 2 3 4 in colors
for i in range(5):
    ax.scatter(subcube_data[i][:,0], 
               subcube_data[i][:,1], 
               subcube_data[i][:,2], marker='o', s=2)
for i in range(5,len(subcube_data)):
    ax.scatter(subcube_data[i][:,0], 
               subcube_data[i][:,1],
                subcube_data[i][:,2],marker='o', s=1, color='black')