From: Waichler, Scott R <Scott.Waichler_at_pnl.gov>

Date: Tue, 10 Mar 2009 11:19:10 -0700

}

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

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.
*