[R] Alternative to interp.surface() offered

From: Waichler, Scott R <Scott.Waichler_at_pnl.gov>
Date: Tue, 10 Mar 2009 11:19:10 -0700


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



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue 10 Mar 2009 - 17:32:36 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 10 Mar 2009 - 18:30:26 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive