[R] Alternative to interp.surface() offered

[Waichler, Scott R]

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