I wanted a simple function for bilinear interpolation on a 2-D grid, and
interp.surface() in the fields package didn't quite suit my needs. In
particular, it requires uniform spacing between grid points. It also
didn't have the "visual" reference frame I was looking for. Here
is an
alternative function, followed by an example.
# A function for bilinear interpolation on a 2-d grid, based on the
# interp.surface() from the fields package and code by Steve Koehler.
# The points of the grid do not have to be uniformly spaced. Looking at
the 2-d
# grid in plan view, the origin is at upper left, so y (row) index
increases
# downward. findInterval() is used to locate the new coordinates on the
grid.
#
# Scott Waichler, scott.waichler at pnl.gov, 03/10/09.
my.interp.surface <- function (obj, loc) {
# obj is a surface object like the list for contour or image.
# loc is a matrix of (x, y) locations
x <- obj$x
y <- 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])
# set weights for out-of-bounds locations to NA
ex[ex < 0 | ex > 1] <- NA
ey[ey < 0 | ey > 1] <- NA
return(
z[cbind(ind.y, ind.x)] * (1 - ex) * (1 - ey) + # upper
left
z[cbind(ind.y + 1, ind.x)] * (1 - ex) * ey + # lower
left
z[cbind(ind.y + 1, ind.x + 1)] * ex * ey + # lower
right
z[cbind(ind.y, ind.x + 1)] * ex * (1 - ey) # upper
right
)
}
## # An example.
## # z matrix, y index increasing downwards
## # 4 5 6 7 8
## # 3 4 5 6 7
## # 2 3 4 5 6
## # 1 2 3 4 5
## z.vec <- c(4,5,6,7,8,3,4,5,6,7,2,3,4,5,6,1,2,3,4,5) # "read in"
the
data for the matrix
## x.mat <- 1:5 # x coordinates of the z values
## y.mat <- seq(100, 400, by=100) # y coordinates of the z values
## obj <- list(x = x.mat, y = y.mat, z = matrix(z.vec, ncol=5, byrow=T))
# grid you want to interpolate on
## x.out <- round(runif(6, min = min(x.mat), max = max(x.mat)), 2) # x
for points you want interpolate to
## y.out <- round(runif(6, min = min(y.mat), max = max(y.mat)), 2) # y
for points you want interpolate to
## loc <- cbind(x.out, y.out)
## z.out <- my.interp.surface(obj, loc)
## cat(file="", "x.out = ", loc[,1], "\n",
"y.out = ", loc[,2], "\n",
"z.out = ", round(z.out, 2), "\n")
Regards,
Scott Waichler
Pacific Northwest National Laboratory
Richland, WA 99352 USA
scott.waichler at pnl.gov
