interp {akima}R Documentation

Gridded Bivariate Interpolation for Irregular Data

Description

If ncp is zero, linear interpolation is used in the triangles bounded by data points. Cubic interpolation is done if partial derivatives are used. If extrap is FALSE, z-values for points outside the convex hull are returned as NA. No extrapolation can be performed if ncp is zero.

The interp function handles duplicate (x,y) points in different ways. As default it will stop with an error message. But it can give duplicate points an unique z value according to the parameter duplicate (mean,median or any other user defined function).

The triangulation scheme used by interp works well if x and y have similar scales but will appear stretched if they have very different scales. The spreads of x and y must be within four orders of magnitude of each other for interp to work.

Usage

interp(x, y, z, xo=<<see below>>, yo=<<see below>>, ncp=0, extrap=F)
interp.old(x, y, z, xo=<<see below>>, yo=<<see below>>, ncp=0, extrap=F)
interp.new(x, y, z, xo=<<see below>>, yo=<<see below>>, ncp=0, extrap=F)

Arguments

x vector of x-coordinates of data points. Missing values are not accepted.
y vector of y-coordinates of data points. Missing values are not accepted.
z vector of z-coordinates of data points. Missing values are not accepted.

x, y, and z must be the same length and may contain no fewer than four points. The points of x and y cannot be collinear, i.e, they cannot fall on the same line (two vectors x and y such that y = ax + b for some a, b will not be accepted). interp is meant for cases in which you have x, y values scattered over a plane and a z value for each. If, instead, you are trying to evaluate a mathematical function, or get a graphical interpretation of relationships that can be described by a polynomial, try outer().

xo vector of x-coordinates of output grid. The default is 40 points evenly spaced over the range of x. If extrapolation is not being used (extrap=F, the default), xo should have a range that is close to or inside of the range of x for the results to be meaningful.
yo vector of y-coordinates of output grid. The default is 40 points evenly spaced over the range of y. If extrapolation is not being used (extrap=F, the default), yo should have a range that is close to or inside of the range of y for the results to be meaningful.
ncp number of additional points to be used in computing partial derivatives at each data point. ncp must be either 0 (partial derivatives are not used), or at least 2 but smaller than the number of data points (and smaller than 25). This option is only supported by interp.old.
extrap logical flag: should extrapolation be used outside of the convex hull determined by the data points?
duplicate indicates how to handle duplicate data points. Possible values are "error" - produces an error message, "strip" - remove duplicate z values, "mean","median","user" - calculate mean , median or user defined function of duplicate z values.
dupfun this function is applied to duplicate points if duplicate="user"

Value

list with 3 components:

x vector of x-coordinates of output grid, the same as the input argument xo, if present. Otherwise, a vector 40 points evenly spaced over the range of the input x.
y vector of y-coordinates of output grid, the same as the input argument yo, if present. Otherwise, a vector 40 points evenly spaced over the range of the input x.
z matrix of fitted z-values. The value z[i,j] is computed at the x,y point x[i], y[j]. z has dimensions length(x) times length(y) (length(xo) times length(yo)).

Note

interp is a wrapper for the two versions interp.old (it uses (almost) the same Fortran code from Akima 1978 as the S-Plus version) and interp.new (it is based on new Fortran code from Akima 1996). For linear interpolation the old version is choosen, but spline interpolation is done by the new version.

At the moment interp.new ignores ncp and does only bicubic spline interpolation.

The resulting structure is suitable for input to the functions contour and image. Check the requirements of these functions when choosing values for xo and yo.

References

Akima, H. (1978). A Method of Bivariate Interpolation and Smooth Surface Fitting for Irregularly Distributed Data Points. ACM Transactions on Mathematical Software, 4, 148-164.

Akima, H. (1996). Algorithm 761: scattered-data surface fitting that has the accuracy of a cubic polynomial. ACM Transactions on Mathematical Software, 22, 362–371

See Also

contour, image, approx, spline, outer, expand.grid.

Examples

data(akima)
# linear interpolation
akima.li <- interp(akima$x, akima$y, akima$z)
image(akima.li$x,akima.li$y,akima.li$z)
contour(akima.li$x,akima.li$y,akima.li$z,add=T)
points(akima$x,akima$y)

# increase smoothness
akima.smooth <- interp(akima$x, akima$y, akima$z,
      xo=seq(0,25, length=100),  yo=seq(0,20, length=100))
image(akima.smooth$x,akima.smooth$y,akima.smooth$z)
contour(akima.smooth$x,akima.smooth$y,akima.smooth$z,add=T)
points(akima$x,akima$y)
# use triangulation library to
# show underlying triangulation:
if(library(tripack, logical.return=T))
  plot(tri.mesh(akima),add=T,lty="dashed")

# use only 15 points (interpolation only within convex hull!)
akima.part <- interp(akima$x[1:15],akima$y[1:15],akima$z[1:15])
image(akima.part$x,akima.part$y,akima.part$z)
contour(akima.part$x,akima.part$y,akima.part$z,add=T)
points(akima$x[1:15],akima$y[1:15])

# spline interpolation, use 5 points to calculate derivatives
# interp gives `linear interpolation not yet implemented with interp.new()'
akima.spl <- interp.old(akima$x, akima$y, akima$z,
      xo=seq(0,25, length=100),  yo=seq(0,20, length=100),ncp=5)
image(akima.spl$x,akima.spl$y,akima.spl$z)
contour(akima.spl$x,akima.spl$y,akima.spl$z,add=T)
points(akima$x,akima$y)

# example with duplicate points
data(airquality)
air <- airquality[(!is.na(airquality$Temp) & 
                   !is.na(airquality$Ozone) & 
                   !is.na(airquality$Solar.R)),]
# gives an error:
air.ip <- interp(air$Temp,air$Solar.R,air$Ozone)
# use mean of duplicate points:
air.ip <- interp(air$Temp,air$Solar.R,air$Ozone,duplicate="mean")

[Package Contents]