CODEWITHSUNDEEP
3dgridscipyinterpolation
Aisha Down
Aisha Down

Reputation: 41

3d Interpolation in Scipy--a density grid

I have a set of coordinates representing the 3d positions of a number of objects. These points are taken from a simulation of a 3-d cube. What I need to do is to make a 3d grid on the cube, and then to assign these points to their correct places on the grid, so that I can then find the density of objects in each section of the grid. I've been searching the interpolation and grid documentation (http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata), for instance, but am unsure what to make of it because I don't have a function associated with these data points. What I was thinking of doing was using an if statement: given an array of points=[x1,y1,z1],[x2,y2,z2],etc, if points[i][0]-gridpoint1[0]<1: if points[i][1]-gridpoint1[1]<1: if points[i][2]-gridpoint1[2]<1, points[i]=bin1[i], where bin1 would be a premade array of zeros. However, I think I'd have to run this for every gridpoint on the grid (the gridpoint would be in the center of every bin), and then to figure out how many non-zero elements there were in each bin, which I am also not sure how to do. I have a feeling I can use some sort of function in scipy to do this entire task more easily, but I am still unsure how to get there. Thanks so much in advance for your help!

Upvotes: 4

Views: 1360

Answers (1)

pv.
pv.

Reputation: 35125

If I understand correctly, you have coordinates (x[i], y[i], z[i]), i = 0, ..., N-1, and you want to count how many of them end up in a given grid cell 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 find out which grid cell each of the points is in, this can be done with 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]

Upvotes: 4

Related Questions