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

[Michael Sumner]

We do a similar thing in 2D using tabulate. It assumes that the coordinates
are known not to fall outside the grid:

## xx are the x-coordinates, similarly yy
## orig is the x/y offset of the grid
## scal is the x/y width and height of grid cells
## x/ydim is the dimensions of the grid

cps <- ceiling(cbind((xx - orig[1])/scl[1], (yy - orig[2])/scl[2]))
tps <- tabulate((cps[, 1] - 1) * ydim + cps[, 2], xdim * ydim) 

## then shape the counts as a matrix, appropriately aligned:
mps <- matrix(tps, ydim, xdim)
z <- t(mps)

I'm sure you could adapt this to the 3D case. That fragment comes from the
function countPoints (called by tripGrid) from the package "trip". 

HTH
Mike

-----Original Message-----
From: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Waichler, Scott R
Sent: Friday, 17 August 2007 8:01 AM
To: R-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Need fast method to find indices of cells in a 3-D
gridthat contain (x, y, z) coordinates

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

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo