[R-sig-Geo] Classify x, y, z points based on upper-boundary surfaces

[Waichler, Scott R]

Hi,

I have surfaces in ascii grid format that represent the elevation of the
upper boundary of material types.  I have multiple ascii grid files, one
for each surface.  The ascii grid files have the same x, y origin and
discretization, but the elevations are continuous variables.  I would
like to create a 3D grid of material type, assigning a type for each x,
y, z location I define.  The correct type depends on where the point
fits in the "sandwich."  So it's a kind of overlay problem.

library(sp)

# Create some ascii grid files
asciigridfn <- c("surf1", "surf2", "surf3")
surfaces <- character(3)
surfaces[1] <- "ncols         4
nrows         3
xllcorner     0
yllcorner     0
cellsize      1
NODATA_value  -999
9.56  9.52  9.11  9.74
9.91 10.00  9.61  9.46
9.13  9.59  9.93  9.55
"
surfaces[2] <- "ncols         4
nrows         3
xllcorner     0
yllcorner     0
cellsize      1
NODATA_value  -999
6.81 6.06 6.42 6.41
6.93 -999 6.84 6.55
6.98 6.15 6.43 6.66
"
surfaces[3] <- "ncols         4
nrows         3
xllcorner     0
yllcorner     0
cellsize      1
NODATA_value  -999
3.82 3.83 3.92 3.28
3.74 3.28 3.26 3.42
3.88 3.69 3.34 3.91
"

# I write temporary files because I couldn't get textConnection() to
work here.
g <- list()
for(i in 1:3) {
  # Write the file
  cat(file=asciigridfn[i], surfaces[i])

  # Read it with read.asciigrid()
  g[[i]] <- read.asciigrid(asciigridfn[i])
}
unlink(asciigridfn)  # delete the files

# Define x, y, z coordinates for points to classify as type 1, 2, 3 or
NA.
# I seek a method to do the classification at each point (results
below).  
# Is there a way with sp?  It's a 3D overlay with some kind of melding
of list g.
# I can only think of a column-by-column analysis with findInterval().  
# My 2-D grids of surfaces have ~1e7 cells.
p1 <- c(0.5, 0.5, 3.5) # type 3
p2 <- c(0.5, 0.5, 3.9) # type 2
p3 <- c(1.5, 1.5, 11)  # type is NA; point located above highest surface
p4 <- c(1.5, 1.5, 4)   # type 1 (no type 2 at this x,y location)
p5 <- c(1.5, 1.5, 3)   # type 3

Thanks,
Scott Waichler
Pacific Northwest National Laboratory
scott.waichler at pnl.gov