[R] Alternative to interp.surface() offered

[Waichler, Scott R]

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@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