akima (version 0.5-7)

interp: Gridded Bivariate Interpolation for Irregular Data

Description

These functions implement bivariate interpolation onto a grid for irregularly spaced input data. Bilinear or bicubic spline interpolation is applied using different versions of algorithms from Akima.

Usage

interp(x, y, z, xo=seq(min(x), max(x), length = 40),
       yo=seq(min(y), max(y), length = 40),
       linear = TRUE, extrap=FALSE, duplicate = "error", dupfun = NULL, ncp = NULL)
interp.old(x, y, z, xo= seq(min(x), max(x), length = 40),
           yo=seq(min(y), max(y), length = 40), ncp = 0,
           extrap=FALSE, duplicate = "error", dupfun = NULL)
interp.new(x, y, z, xo = seq(min(x), max(x), length = 40),
           yo = seq(min(y), max(y), length = 40), linear = FALSE,
            ncp = NULL, extrap=FALSE, duplicate = "error", dupfun = NULL)

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

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=FALSE, the default), xo should have a range that is close t

yo

vector of y-coordinates of output grid; analogous to xo, see above.

linear

logical -- indicating wether linear or spline interpolation should be used. supersedes old ncp parameter

ncp

deprecated, use parameter linear. Now only used by interp.old().

meaning was: number of additional points to be used in computing partial derivatives at each data point. ncp must be either

extrap

logical flag: should extrapolation be used outside of the convex hull determined by the data points?

duplicate

character string indicating how to handle duplicate data points. Possible values are [object Object],[object Object],[object Object]

dupfun

a function, applied to duplicate points if duplicate= "user".

Value

• list with 3 components:
• x,y
  vectors of x- and y- coordinates of output grid, the same as the input argument xo, or yo, if present. Otherwise, their default, 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 xo[i], yo[j]. z has dimensions length(xo) times length(yo).

Details

If linear is TRUE (default), linear interpolation is used in the triangles bounded by data points. Cubic interpolation is done if linear is set to FALSE. If extrap is FALSE, z-values for points outside the convex hull are returned as NA. No extrapolation can be performed for the linear case.

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.

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, aspline, outer, expand.grid.

Examples

data(akima)
plot(y ~ x, data = akima, main = "akima example data")
with(akima, text(x, y, formatC(z,dig=2), adj = -0.1))

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

## increase smoothness (using finer grid):
akima.smooth <-
 with(akima, interp(x, y, z, xo=seq(0,25, length=100),
                    yo=seq(0,20, length=100)))
image  (akima.smooth, main = "interp(<akima data>, *) on finer grid")
contour(akima.smooth, add = TRUE, col = "thistle")
points(akima, pch = 3, cex = 2, col = "blue")
# use triangulation package to show underlying triangulation:
if(library(tripack, logical.return=TRUE))
   plot(tri.mesh(akima), add=TRUE, lty="dashed")
# use only 15 points (interpolation only within convex hull!)
akima.part <- with(akima, interp(x[1:15], y[1:15], z[1:15]))
image(akima.part)
title("interp() on subset of only 15 points")
contour(akima.part, add=TRUE)
points(akima$x[1:15],akima$y[1:15], col = "blue")

## spline interpolation, two variants (AMS 526 "Old", AMS 761 "New")
## -----------------------------------------------------------------
## "Old": use 5 points to calculate derivatives -> many NAs
akima.sO <- interp.old(akima$x, akima$y, akima$z,
           xo=seq(0,25, length=100),  yo=seq(0,20, length=100), ncp=5)
table(is.na(akima.sO$z)) ## 3990 NA's; = 40 \%
akima.sO <- with(akima,
   interp.old(x,y,z, xo=seq(0,25, length=100), yo=seq(0,20, len=100), ncp = 4))
sum(is.na(akima.sO$z)) ## still 3429
image  (akima.sO, main = "interp.old(*, ncp = 4) [almost useless]")
contour(akima.sO, add = TRUE)

## "New:"
akima.spl <- with(akima, interp.new(x,y,z, xo=seq(0,25, length=100),
                                           yo=seq(0,20, length=100)))
## equivalent call via setting linear=FALSE in interp():
akima.spl <- with(akima, interp(x,y,z, xo=seq(0,25, length=100),
                                       yo=seq(0,20, length=100),
                                linear=FALSE))

contour(akima.spl, main = "smooth  interp(*, linear = FALSE)")
points(akima)

full.pal <- function(n) hcl(h = seq(340, 20, length = n))
cool.pal <- function(n) hcl(h = seq(120, 0, length = n) + 150)
warm.pal <- function(n) hcl(h = seq(120, 0, length = n) - 30)

filled.contour(akima.spl, color.palette = full.pal,
        plot.axes = { axis(1); axis(2);
                      title("smooth  interp(*, linear = FALSE)");
                      points(akima, pch = 3, col= hcl(c=100, l = 20))})
# no extrapolation!

## example with duplicate points :

data(airquality)
air <- subset(airquality,
              !is.na(Temp) & !is.na(Ozone) & !is.na(Solar.R))
# gives an error {duplicate ..}:
try( air.ip <- interp(air$Temp,air$Solar.R,air$Ozone, linear=FALSE) )
# use mean of duplicate points:
air.ip <- with(air, interp(Temp, Solar.R, log(Ozone), duplicate = "mean",
                           linear = FALSE))
image(air.ip, main = "Airquality: Ozone vs. Temp and Solar.R")
with(air, points(Temp, Solar.R))