3d Interpolation in Scipy - density grid

Question

I have a set of coordinates, which is a three-dimensional position of a number of objects. These points are taken from a 3D cube simulation. I need to make a 3D grid on a cube, and then assign these points to their correct places in the grid, in order to then find the density of objects in each section of the grid. For example, I was looking for interpolation and grid documentation (http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata), but I'm not sure what to do, because I do not have a function associated with these data points. What I was going to do is use the if statement: an array of points = [x1, y1, z1], [x2, y2, z2], etc., if the points are [i] [0] -gridpoint1 [0] & lt 1: if points [i] [1] -gridpoint1 [1] <1: if points [i] [2] -gridpoint1 [2] <1, points [i] = bin1 [i] where bin1 will be ready an array of zeros. However, I think that I will have to run this for each grid point on the grid (the grid point will be in the center of each bin), and then find out how many non-zero elements were in each box, which I'm also not sure how to do . I have a feeling that I can use some kind of function in scipy to make this task more easy, but I'm still not sure how to get there. Many thanks for your help!

Answer 1

If I understand correctly, you have the coordinates (x [i], y [i], z [i]), i = 0, ..., N-1, and you want to calculate how many of them end in the given cell grids in a 3D cube?

This can be done using numpy.histogramdd :

 import numpy as np # some random points in a cube x = np.random.rand(100) y = np.random.rand(100) z = np.random.rand(100); # specify grid in a cube xgrid = np.linspace(0.0, 1.0, 5) ygrid = np.linspace(0.0, 1.0, 6) zgrid = np.linspace(0.0, 1.0, 7) # compute counts counts, edges = np.histogramdd(np.c_[x, y, z], bins=(xgrid, ygrid, zgrid)) print counts.shape # -> (4, 5, 6) # number of points in the box # xgrid[1] <= x < xgrid[2] && ygrid[2] <= y < ygrid[3] && zgrid[0] <= z < zgrid[1] print counts[1, 2, 0] 

If you want to know which mesh cell is at each point, this can be done using searchsorted :

 ix = np.searchsorted(xgrid, x) - 1 iy = np.searchsorted(ygrid, y) - 1 iz = np.searchsorted(zgrid, z) - 1 # point (x[3], y[3], z[3]) is in the following grid cell print ix[3], iy[3], iz[3]