This is a correction to a post from 3/10/09. I just wanted to get this in the
archive. It is the same thread as
http://www.mail-archive.com/r-help at r-project.org/msg48762.html
Thanks to Matt Oliver for bringing this to my attention.
The correct code for my.interp.surface() follows.
# A function for bilinear interpolation on a 2-d grid, based on
# interp.surface() from the fields package and code by Steve Koehler.
# The points of the grid do not have to be uniformly spaced.
# findInterval() is used to locate on the grid.
my.interp.surface <- function (obj, loc) {
# obj is a list with known x, y, z.
# loc a matrix of new x, y locations at which you want z values.
x <- sort(unique(obj$x))
y <- sort(unique(obj$y))
x.new <- loc[,1]
y.new <- loc[,2]
z <- obj$z
ind.x <- findInterval(x.new, x, all.inside=T)
ind.y <- findInterval(y.new, y, all.inside=T)
ex <- (x.new - x[ind.x]) / (x[ind.x + 1] - x[ind.x])
ey <- (y.new - y[ind.y]) / (y[ind.y + 1] - y[ind.y])
ex[ex < 0 | ex > 1] <- NA
ey[ey < 0 | ey > 1] <- NA
return(
z[cbind(ind.x, ind.y) ] * (1 - ex) * (1 - ey) + # upper left
z[cbind(ind.x, ind.y + 1) ] * (1 - ex) * ey + # lower left
z[cbind(ind.x + 1, ind.y + 1)] * ex * ey + # lower right
z[cbind(ind.x + 1, ind.y) ] * ex * (1 - ey) # upper right
)
}
Scott Waichler
Pacific Northwest National Laboratory
scott.waichler at pnnl.gov