[R-sig-Geo] Need fast method to find indices of cells in a 3-D grid that contain (x, y, z) coordinates

[Waichler, Scott R]

I am looking for a fast method that will find the indices of cells of a
3-D grid that contain a set of (x,y,z) coordinates.  The grid has I,J,K
number of cells in the x-, y-, and z- directions, respectively.  The
cells are numbered with integers 1:(I*J*K).  I have multiple sets of
coordinates (points in 3-D space) for which I want to rapidly look up
the numbers of the cells that contain them.  Here is an example of what
I am trying to do, with an explicit loop that evaluates each point
separately.  Is there a faster way that can do the whole set of points
at once?

# Define a 10x10x5 3-D grid in terms of lower and upper bounds of cells
I <- 10  # number of cells in x-direction
J <- 10  # number of cells in y-direction
K <- 5   # number of cells in z-direction
x.lower.bounds <- 0:(I-1)
x.upper.bounds <- 1:I
y.lower.bounds <- 0:(J-1)
y.upper.bounds <- 1:J
z.lower.bounds <- 0:(K-1)
z.upper.bounds <- 1:K
num.grid.cells <- I*J*K

# Now repeat the above vectors to create new vectors for all cells in
the 3-D grid.
# I know this could be faster with vectorization, but I don't need speed
here :).
x.lb <- x.ub <- y.lb <- y.ub <- z.lb <- z.ub <- numeric(num.grid.cells)
n <- 0
for(k in 1:K) {
  for(j in 1:J) {
    for(i in 1:I) {
      n <- n + 1
      x.lb[n] <- x.lower.bounds[i]
      x.ub[n] <- x.upper.bounds[i]
      y.lb[n] <- y.lower.bounds[j]
      y.ub[n] <- y.upper.bounds[j]
      z.lb[n] <- z.lower.bounds[k]
      z.ub[n] <- z.upper.bounds[k]
    }
  }
}
  
# Define some points for which I want to look up the grid cell numbers:
num.points <- 3
x <- runif(num.points, 0, I)
y <- runif(num.points, 0, J)
z <- runif(num.points, 0, K)

# Loop thru the points and find the grid cell number of each.  Is there
a 
# way to vectorize this and do all points at once, either explicitly
with base 
# R or in some other way with the help of a package?
for(n in 1:num.points) {
 grid.cell.number <- which(x.lb <= x[n] & x.ub >= x[n] &
                           y.lb <= y[n] & y.ub >= y[n] &
                           z.lb <= z[n] & z.ub >= z[n])
 cat("The grid cell number for point", x[n], y[n], z[n], "is",
grid.cell.number, "\n")
}

Regards,
Scott Waichler
Pacific Northwest National Laboratory
Richland, Washington, USA
scott.waichler _at_ pnl.gov