| Title: | Interpolation Methods |
|---|---|
| Description: | Bivariate data interpolation on regular and irregular grids, either linear or using splines are the main part of this package. It is intended to provide FOSS replacement functions for the ACM licensed akima::interp and tripack::tri.mesh functions. Linear interpolation is implemented in interp::interp(..., method="linear"), this corresponds to the call akima::interp(..., linear=TRUE) which is the default setting and covers most of akima::interp use cases in depending packages. A re-implementation of Akimas irregular grid spline interpolation (akima::interp(..., linear=FALSE)) is now also available via interp::interp(..., method="akima"). Estimators for partial derivatives are now also available in interp::locpoly(), these are a prerequisite for the spline interpolation. The basic part is a GPLed triangulation algorithm (sweep hull algorithm by David Sinclair) providing the starting point for the irregular grid interpolator. As side effect this algorithm is also used to provide replacements for almost all functions of the tripack package which also suffers from the same ACM license restrictions. All functions are designed to be backward compatible with their akima / tripack counterparts. |
| Authors: | Albrecht Gebhardt [aut, cre, cph], Roger Bivand [aut], David Sinclair [aut, cph] (author of the shull library) |
| Maintainer: | Albrecht Gebhardt <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.1-6 |
| Built: | 2024-11-22 06:59:10 UTC |
| Source: | CRAN |
Interpolation of values given regular or irregular gridded
data sets containing coordinates and function values
is (will be) available through this package. As this
interpolation is (for the irregular gridded data case) based on
trianglation of the data locations also triangulation functions are
implemented. Moreover the (not yet finished) spline interpolation
needs estimators for partial derivates, these are also made available
to the end user for direct use.
The interpolation use can be divided by the used method into piecewise linear (finished in 1_0.27) and spline (not yet finished) interpolation and by input and output settings into gridded and pointwise setups.
This package is a FOSS replacement for the ACM licensed packages
akima and tripack. The function calls are backward
compatible.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
Maintainer: Albrecht Gebhardt <[email protected]>
interp, tri.mesh,
voronoi.mosaic, locpoly
akima is a list with components x, y and z which
represents a smooth surface of z values at selected points
irregularly distributed in the x-y plane.
The data was taken from a study of waveform distortion in electronic circuits, described in: Hiroshi Akima, "A Method of Bivariate Interpolation and Smooth Surface Fitting Based on Local Procedures", CACM, Vol. 17, No. 1, January 1974, pp. 18-20.
Hiroshi Akima, "A Method of Bivariate Interpolation and Smooth Surface Fitting for Irregularly Distributed Data Points", ACM Transactions on Mathematical Software, Vol. 4, No. 2, June 1978, pp. 148-159. Copyright 1978, Association for Computing Machinery, Inc., reprinted by permission.
## Not run: library(rgl) data(akima) # data rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red") rgl.bbox() # bivariate linear interpolation # interp: akima.li <- interp(akima$x, akima$y, akima$z, xo=seq(min(akima$x), max(akima$x), length = 100), yo=seq(min(akima$y), max(akima$y), length = 100)) # interp surface: rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5)) # interpp: akima.p <- interpp(akima$x, akima$y, akima$z, runif(200,min(akima$x),max(akima$x)), runif(200,min(akima$y),max(akima$y))) # interpp points: rgl.points(akima.p$x,akima.p$z , akima.p$y,size=4,color="yellow") # bivariate spline interpolation # data rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red") rgl.bbox() # bivariate cubic spline interpolation # interp: akima.si <- interp(akima$x, akima$y, akima$z, xo=seq(min(akima$x), max(akima$x), length = 100), yo=seq(min(akima$y), max(akima$y), length = 100), linear = FALSE, extrap = TRUE) # interp surface: rgl.surface(akima.si$x,akima.si$y,akima.si$z,color="green",alpha=c(0.5)) # interpp: akima.sp <- interpp(akima$x, akima$y, akima$z, runif(200,min(akima$x),max(akima$x)), runif(200,min(akima$y),max(akima$y)), linear = FALSE, extrap = TRUE) # interpp points: rgl.points(akima.sp$x,akima.sp$z , akima.sp$y,size=4,color="yellow") ## End(Not run)## Not run: library(rgl) data(akima) # data rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red") rgl.bbox() # bivariate linear interpolation # interp: akima.li <- interp(akima$x, akima$y, akima$z, xo=seq(min(akima$x), max(akima$x), length = 100), yo=seq(min(akima$y), max(akima$y), length = 100)) # interp surface: rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5)) # interpp: akima.p <- interpp(akima$x, akima$y, akima$z, runif(200,min(akima$x),max(akima$x)), runif(200,min(akima$y),max(akima$y))) # interpp points: rgl.points(akima.p$x,akima.p$z , akima.p$y,size=4,color="yellow") # bivariate spline interpolation # data rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red") rgl.bbox() # bivariate cubic spline interpolation # interp: akima.si <- interp(akima$x, akima$y, akima$z, xo=seq(min(akima$x), max(akima$x), length = 100), yo=seq(min(akima$y), max(akima$y), length = 100), linear = FALSE, extrap = TRUE) # interp surface: rgl.surface(akima.si$x,akima.si$y,akima.si$z,color="green",alpha=c(0.5)) # interpp: akima.sp <- interpp(akima$x, akima$y, akima$z, runif(200,min(akima$x),max(akima$x)), runif(200,min(akima$y),max(akima$y)), linear = FALSE, extrap = TRUE) # interpp points: rgl.points(akima.sp$x,akima.sp$z , akima.sp$y,size=4,color="yellow") ## End(Not run)
akima474 is a list with vector components x, y and a matrix z which
represents a smooth surface of z values at the points
of a regular grid spanned by the vectors x and y.
Hiroshi Akima, Bivariate Interpolation and Smooth Surface Fitting Based on Local Procedures [E2], Communications of ACM, Vol. 17, No. 1, January 1974, pp. 26-30
## Not run: library(rgl) data(akima474) # data rgl.spheres(akima474$x,akima474$z , akima474$y,0.5,color="red") rgl.bbox() # bivariate linear interpolation # interp: akima474.li <- interp(akima474$x, akima474$y, akima474$z, xo=seq(min(akima474$x), max(akima474$x), length = 100), yo=seq(min(akima474$y), max(akima474$y), length = 100)) # interp surface: rgl.surface(akima474.li$x,akima474.li$y,akima474.li$z,color="green",alpha=c(0.5)) # interpp: akima474.p <- interpp(akima474$x, akima474$y, akima474$z, runif(200,min(akima474$x),max(akima474$x)), runif(200,min(akima474$y),max(akima474$y))) # interpp points: rgl.points(akima474.p$x,akima474.p$z , akima474.p$y,size=4,color="yellow") # bivariate spline interpolation # data rgl.spheres(akima474$x,akima474$z , akima474$y,0.5,color="red") rgl.bbox() # bivariate cubic spline interpolation # interp: akima474.si <- interp(akima474$x, akima474$y, akima474$z, xo=seq(min(akima474$x), max(akima474$x), length = 100), yo=seq(min(akima474$y), max(akima474$y), length = 100), linear = FALSE, extrap = TRUE) # interp surface: rgl.surface(akima474.si$x,akima474.si$y,akima474.si$z,color="green",alpha=c(0.5)) # interpp: akima474.sp <- interpp(akima474$x, akima474$y, akima474$z, runif(200,min(akima474$x),max(akima474$x)), runif(200,min(akima474$y),max(akima474$y)), linear = FALSE, extrap = TRUE) # interpp points: rgl.points(akima474.sp$x,akima474.sp$z , akima474.sp$y,size=4,color="yellow") ## End(Not run)## Not run: library(rgl) data(akima474) # data rgl.spheres(akima474$x,akima474$z , akima474$y,0.5,color="red") rgl.bbox() # bivariate linear interpolation # interp: akima474.li <- interp(akima474$x, akima474$y, akima474$z, xo=seq(min(akima474$x), max(akima474$x), length = 100), yo=seq(min(akima474$y), max(akima474$y), length = 100)) # interp surface: rgl.surface(akima474.li$x,akima474.li$y,akima474.li$z,color="green",alpha=c(0.5)) # interpp: akima474.p <- interpp(akima474$x, akima474$y, akima474$z, runif(200,min(akima474$x),max(akima474$x)), runif(200,min(akima474$y),max(akima474$y))) # interpp points: rgl.points(akima474.p$x,akima474.p$z , akima474.p$y,size=4,color="yellow") # bivariate spline interpolation # data rgl.spheres(akima474$x,akima474$z , akima474$y,0.5,color="red") rgl.bbox() # bivariate cubic spline interpolation # interp: akima474.si <- interp(akima474$x, akima474$y, akima474$z, xo=seq(min(akima474$x), max(akima474$x), length = 100), yo=seq(min(akima474$y), max(akima474$y), length = 100), linear = FALSE, extrap = TRUE) # interp surface: rgl.surface(akima474.si$x,akima474.si$y,akima474.si$z,color="green",alpha=c(0.5)) # interpp: akima474.sp <- interpp(akima474$x, akima474$y, akima474$z, runif(200,min(akima474$x),max(akima474$x)), runif(200,min(akima474$y),max(akima474$y)), linear = FALSE, extrap = TRUE) # interpp points: rgl.points(akima474.sp$x,akima474.sp$z , akima474.sp$y,size=4,color="yellow") ## End(Not run)
This function extracts a list of arcs from a triangulation object
created by tri.mesh.
arcs(tri.obj)arcs(tri.obj)
tri.obj |
object of class |
This function acesses the arcs component of a triangulation
object returned by tri.mesh and extracts the arcs
contained in this triangulation. This is e.g. used for plotting.
A matrix with two columns "from" and "to" containing the
indices of points connected by the arc with the corresponding row index.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
data(franke) tr <- tri.mesh(franke$ds3) arcs(tr)data(franke) tr <- tri.mesh(franke$ds3) arcs(tr)
This function returns a list containing the areas of each triangle of
a triangulation object created by tri.mesh.
area(tri.obj)area(tri.obj)
tri.obj |
object of class |
This function acesses the cclist component of a triangulation
object returned by tri.mesh and extracts the areas
of the triangles contained in this triangulation.
A vector containing the area values.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
data(franke) tr <- tri.mesh(franke$ds3) area(tr)data(franke) tr <- tri.mesh(franke$ds3) area(tr)
The function returns a list of points which smoothly interpolate given data points, similar to a curve drawn by hand.
aspline(x, y = NULL, xout, n = 50, ties = mean, method = "improved", degree = 3) aSpline(x, y, xout, method = "improved", degree = 3)aspline(x, y = NULL, xout, n = 50, ties = mean, method = "improved", degree = 3) aSpline(x, y, xout, method = "improved", degree = 3)
x, y
|
vectors giving the coordinates of the points to be
interpolated. Alternatively a single plotting structure can be
specified: see |
xout |
an optional set of values specifying where interpolation is to take place. |
n |
If |
ties |
Handling of tied |
method |
either |
degree |
if improved algorithm is selected: degree of the polynomials for the interpolating function |
The original algorithm is based on a piecewise function composed of a set of polynomials, each of degree three, at most, and applicable to successive interval of the given points. In this method, the slope of the curve is determined at each given point locally by fitting a third degree polynomial to four consecutive points. Each polynomial representing a portion of the curve between a pair of given points is determined by the coordinates of and the slopes at the points. The data set is prolonged below and above minimum and maximum x values to enable estimation of derivatives at the boundary. The improved algorithm uses polynomials of degree two and one at the boundary. Additionally four overlapping sequences of points are used for the estimation via a residual based weighting scheme.
x |
x coordinates of the interpolated data as given by 'xout' or 'n'. |
y |
interpolated y values. |
'aspline' is a wrapper call for the underlying Rcpp function 'aSpline' which could also be called directly with 'x' and 'y' arguments if 'xout' is given and no 'ties' argument is needed.
This is a reimplementation of Akimas algorithms (original and improved version). It is only based on the original articles. It does not involve or resemble the Fortran code associated with those articles. For this reason results may differ slightly because different expressions can result in different numerical errors.
This code is under GPL in contrast to original Fortran code as provided in package 'akima'.
The function arguments are identical to the call in package 'akima', only the 'method' argument has its default now set to 'improved'.
Albrecht Gebhardt <[email protected]>, Thomas Petzold <[email protected]>
Akima, H. (1970) A new method of interpolation and smooth curve fitting based on local procedures, J. ACM 17(4), 589-602
Akima, H. (1991) A Method of Univariate Interpolation that Has the Accuracy of a Third-degree Polynomial. ACM Transactions on Mathematical Software, 17(3), 341-366.
## regular spaced data x <- 1:10 y <- c(rnorm(5), c(1,1,1,1,3)) xnew <- seq(-1, 11, 0.1) plot(x, y, ylim=c(-3, 3), xlim=range(xnew)) ## stats::spline() for comparison lines(spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200), col="blue") lines(aspline(x, y, xnew, method="original"), col="red") lines(aspline(x, y, xnew, method="improved"), col="black", lty="dotted") lines(aspline(x, y, xnew, method="improved", degree=10), col="green", lty="dashed") ## irregular spaced data x <- sort(runif(10, max=10)) y <- c(rnorm(5), c(1,1,1,1,3)) xnew <- seq(-1, 11, 0.1) plot(x, y, ylim=c(-3, 3), xlim=range(xnew)) ## stats::spline() for comparison lines(spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200), col="blue") lines(aspline(x, y, xnew, method="original"), col="red") lines(aspline(x, y, xnew, method="improved"), col="black", lty="dotted") lines(aspline(x, y, xnew, method="improved", degree=10), col="green", lty="dashed") ## an example of Akima, 1991 x <- c(-3, -2, -1, 0, 1, 2, 2.5, 3) y <- c( 0, 0, 0, 0, -1, -1, 0, 2) plot(x, y, ylim=c(-3, 3)) ## stats::spline() for comparison lines(spline(x, y, n=200), col="blue") lines(aspline(x, y, n=200, method="original"), col="red") lines(aspline(x, y, n=200, method="improved"), col="black", lty="dotted") lines(aspline(x, y, n=200, method="improved", degree=10), col="green", lty="dashed")## regular spaced data x <- 1:10 y <- c(rnorm(5), c(1,1,1,1,3)) xnew <- seq(-1, 11, 0.1) plot(x, y, ylim=c(-3, 3), xlim=range(xnew)) ## stats::spline() for comparison lines(spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200), col="blue") lines(aspline(x, y, xnew, method="original"), col="red") lines(aspline(x, y, xnew, method="improved"), col="black", lty="dotted") lines(aspline(x, y, xnew, method="improved", degree=10), col="green", lty="dashed") ## irregular spaced data x <- sort(runif(10, max=10)) y <- c(rnorm(5), c(1,1,1,1,3)) xnew <- seq(-1, 11, 0.1) plot(x, y, ylim=c(-3, 3), xlim=range(xnew)) ## stats::spline() for comparison lines(spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200), col="blue") lines(aspline(x, y, xnew, method="original"), col="red") lines(aspline(x, y, xnew, method="improved"), col="black", lty="dotted") lines(aspline(x, y, xnew, method="improved", degree=10), col="green", lty="dashed") ## an example of Akima, 1991 x <- c(-3, -2, -1, 0, 1, 2, 2.5, 3) y <- c( 0, 0, 0, 0, -1, -1, 0, 2) plot(x, y, ylim=c(-3, 3)) ## stats::spline() for comparison lines(spline(x, y, n=200), col="blue") lines(aspline(x, y, n=200, method="original"), col="red") lines(aspline(x, y, n=200, method="improved"), col="black", lty="dotted") lines(aspline(x, y, n=200, method="improved", degree=10), col="green", lty="dashed")
This is a placeholder function for backward compatibility with packaga akima.
In its current state it simply calls the reimplemented Akima algorithm for irregular grids applied to the regular gridded data given.
Later a reimplementation of the original algorithm for regular grids may follow.
bicubic(x, y, z, x0, y0)bicubic(x, y, z, x0, y0)
x |
a vector containing the |
y |
a vector containing the |
z |
a matrix containing the |
x0 |
vector of |
y0 |
vector of |
This function is a call wrapper for backward compatibility with package akima.
Currently it applies Akimas irregular grid splines to regular grids, later a FOSS reimplementation of his regular grid splines may replace this wrapper.
This function produces a list of interpolated points:
x |
vector of |
y |
vector of |
z |
vector of interpolated data |
If you need an output grid, see bicubic.grid.
Use interp for the general case of irregular gridded data!
Akima, H. (1996) Rectangular-Grid-Data Surface Fitting that Has the Accuracy of a Bicubic Polynomial, J. ACM 22(3), 357-361
data(akima474) # interpolate at the diagonal of the grid [0,8]x[0,10] akima.bic <- bicubic(akima474$x,akima474$y,akima474$z, seq(0,8,length=50), seq(0,10,length=50)) plot(sqrt(akima.bic$x^2+akima.bic$y^2), akima.bic$z, type="l")data(akima474) # interpolate at the diagonal of the grid [0,8]x[0,10] akima.bic <- bicubic(akima474$x,akima474$y,akima474$z, seq(0,8,length=50), seq(0,10,length=50)) plot(sqrt(akima.bic$x^2+akima.bic$y^2), akima.bic$z, type="l")
This is a placeholder function for backward compatibility with packaga akima.
In its current state it simply calls the reimplemented Akima algorithm for irregular grids applied to the regular gridded data given.
Later a reimplementation of the original algorithm for regular grids may follow.
bicubic.grid(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)), nx=40,ny=40,dx=NULL,dy=NULL)bicubic.grid(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)), nx=40,ny=40,dx=NULL,dy=NULL)
x |
a vector containing the |
y |
a vector containing the |
z |
a matrix containing the |
xlim |
vector of length 2 giving lower and upper limit for range |
ylim |
vector of length 2 giving lower and upper limit for range of |
nx |
output grid dimension in |
ny |
output grid dimension in |
dx |
output grid spacing in |
dy |
output grid spacing in |
This function is a call wrapper for backward compatibility with package akima.
Currently it applies Akimas irregular grid splines to regular grids, later a FOSS reimplementation of his regular grid splines may replace this wrapper.
This function produces a grid of interpolated points, feasible to be
used directly with image and contour:
x |
vector of |
y |
vector of |
z |
matrix of interpolated data for the output grid. |
Use interp for the general case of irregular gridded data!
Akima, H. (1996) Rectangular-Grid-Data Surface Fitting that Has the Accuracy of a Bicubic Polynomial, J. ACM 22(3), 357-361
data(akima474) # interpolate at a grid [0,8]x[0,10] akima.bic <- bicubic.grid(akima474$x,akima474$y,akima474$z) zmin <- min(akima.bic$z, na.rm=TRUE) zmax <- max(akima.bic$z, na.rm=TRUE) breaks <- pretty(c(zmin,zmax),10) colors <- heat.colors(length(breaks)-1) image(akima.bic, breaks=breaks, col=colors) contour(akima.bic, levels=breaks, add=TRUE)data(akima474) # interpolate at a grid [0,8]x[0,10] akima.bic <- bicubic.grid(akima474$x,akima474$y,akima474$z) zmin <- min(akima.bic$z, na.rm=TRUE) zmax <- max(akima.bic$z, na.rm=TRUE) breaks <- pretty(c(zmin,zmax),10) colors <- heat.colors(length(breaks)-1) image(akima.bic, breaks=breaks, col=colors) contour(akima.bic, levels=breaks, add=TRUE)
This is an implementation of a bilinear interpolating function.
For a point (x0,y0) contained in a rectangle (x1,y1),(x2,y1), (x2,y2),(x1,y2) and x1<x2, y1<y2, the first step is to get z() at locations (x0,y1) and (x0,y2) as convex linear combinations z(x0,y*)=a*z(x1,y*)+(1-a)*z(x2,y*) where a=(x2-x1)/(x0-x1) for y*=y1,y2. In a second step z(x0,y0) is calculated as convex linear combination between z(x0,y1) and z(x0,y2) as z(x0,y1)=b*z(x0,y1)+(1-b)*z(x0,y2) where b=(y2-y1)/(y0-y1).
Finally, z(x0,y0) is a convex linear combination of the z values at the corners of the containing rectangle with weights according to the distance from (x0,y0) to these corners.
The grid lines can be unevenly spaced.
bilinear(x, y, z, x0, y0) BiLinear(x, y, z, x0, y0)bilinear(x, y, z, x0, y0) BiLinear(x, y, z, x0, y0)
x |
a vector containing the |
y |
a vector containing the |
z |
a matrix containing the |
x0 |
vector of |
y0 |
vector of |
This function produces a list of interpolated points:
x |
vector of |
y |
vector of |
z |
vector of interpolated data |
If you need an output grid, see bilinear.grid.
This Fortran function was part of the akima package but not related to any of Akimas algorithms and under GPL. So it could be transfered into the interp package without changes.
BiLinear is a C++ reimplementation, maybe it will replace the Fortran
implementation later, so
its name may change in futrure versions.
Use interpp for the general case of irregular gridded data!
Pascal Getreuer, Linear Methods for Image Interpolation, Image Processing On Line, 2011, http://www.ipol.im/pub/art/2011/g_lmii/article.pdf
data(akima474) # interpolate at the diagonal of the grid [0,8]x[0,10] akima.bil <- bilinear(akima474$x,akima474$y,akima474$z, seq(0,8,length=50), seq(0,10,length=50)) plot(sqrt(akima.bil$x^2+akima.bil$y^2), akima.bil$z, type="l")data(akima474) # interpolate at the diagonal of the grid [0,8]x[0,10] akima.bil <- bilinear(akima474$x,akima474$y,akima474$z, seq(0,8,length=50), seq(0,10,length=50)) plot(sqrt(akima.bil$x^2+akima.bil$y^2), akima.bil$z, type="l")
This is an implementation of a bilinear interpolating function.
For a point (x0,y0) contained in a rectangle (x1,y1),(x2,y1), (x2,y2),(x1,y2) and x1<x2, y1<y2, the first step is to get z() at locations (x0,y1) and (x0,y2) as convex linear combinations z(x0,y*)=a*z(x1,y*)+(1-a)*z(x2,y*) where a=(x2-x1)/(x0-x1) for y*=y1,y2. In a second step z(x0,y0) is calculated as convex linear combination between z(x0,y1) and z(x0,y2) as z(x0,y1)=b*z(x0,y1)+(1-b)*z(x0,y2) where b=(y2-y1)/(y0-y1).
Finally, z(x0,y0) is a convex linear combination of the z values at the corners of the containing rectangle with weights according to the distance from (x0,y0) to these corners.
The grid lines can be unevenly spaced.
bilinear.grid(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)), nx=40,ny=40,dx=NULL,dy=NULL) BiLinear.grid(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)), nx=40,ny=40,dx=NULL,dy=NULL)bilinear.grid(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)), nx=40,ny=40,dx=NULL,dy=NULL) BiLinear.grid(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)), nx=40,ny=40,dx=NULL,dy=NULL)
x |
a vector containing the |
y |
a vector containing the |
z |
a matrix containing the |
xlim |
vector of length 2 giving lower and upper limit for range |
ylim |
vector of length 2 giving lower and upper limit for range of |
nx |
output grid dimension in |
ny |
output grid dimension in |
dx |
output grid spacing in |
dy |
output grid spacing in |
This function produces a grid of interpolated points, feasible to be
used directly with image and contour:
x |
vector of |
y |
vector of |
z |
matrix of interpolated data for the output grid. |
This Fortran function was part of the akima package but not related to any of Akimas algorithms and under GPL. So it could be transfered into the interp package without changes.
BiLinear.grid is a C++ reimplementation, maybe this will replace the
Fortran implementation later. So its name may change in future versions, dont
rely on it currently.
Pascal Getreuer, Linear Methods for Image Interpolation, Image Processing On Line, 2011, http://www.ipol.im/pub/art/2011/g_lmii/article.pdf
data(akima474) # interpolate at a grid [0,8]x[0,10] akima.bil <- bilinear.grid(akima474$x,akima474$y,akima474$z) zmin <- min(akima.bil$z, na.rm=TRUE) zmax <- max(akima.bil$z, na.rm=TRUE) breaks <- pretty(c(zmin,zmax),10) colors <- heat.colors(length(breaks)-1) image(akima.bil, breaks=breaks, col=colors) contour(akima.bil, levels=breaks, add=TRUE)data(akima474) # interpolate at a grid [0,8]x[0,10] akima.bil <- bilinear.grid(akima474$x,akima474$y,akima474$z) zmin <- min(akima.bil$z, na.rm=TRUE) zmax <- max(akima.bil$z, na.rm=TRUE) breaks <- pretty(c(zmin,zmax),10) colors <- heat.colors(length(breaks)-1) image(akima.bil, breaks=breaks, col=colors) contour(akima.bil, levels=breaks, add=TRUE)
This function returns some info about the cells of a voronoi mosaic, including the coordinates of the vertices and the cell area.
cells(voronoi.obj)cells(voronoi.obj)
voronoi.obj |
object of class |
The function calculates the neighbourhood relations between the underlying triangulation and translates it into the neighbourhood relations between the voronoi cells.
retruns a list of lists, one entry for each voronoi cell which contains
cell |
cell index |
center |
cell 'center' |
neighbours |
neighbour cell indices |
nodes |
2 times |
area |
cell area |
outer cells have area=NA, currently also nodes=NA
which is not really useful – to be done later
A. Gebhardt
data(tritest) tritest.vm <- voronoi.mosaic(tritest$x,tritest$y) tritest.cells <- cells(tritest.vm) # higlight cell 12: plot(tritest.vm) polygon(t(tritest.cells[[12]]$nodes),col="green") # put cell area into cell center: text(tritest.cells[[12]]$center[1], tritest.cells[[12]]$center[2], tritest.cells[[12]]$area)data(tritest) tritest.vm <- voronoi.mosaic(tritest$x,tritest$y) tritest.cells <- cells(tritest.vm) # higlight cell 12: plot(tritest.vm) polygon(t(tritest.cells[[12]]$nodes),col="green") # put cell area into cell center: text(tritest.cells[[12]]$center[1], tritest.cells[[12]]$center[2], tritest.cells[[12]]$area)
This function plots circles at given locations with given radii.
circles(x, y, r, ...)circles(x, y, r, ...)
x |
vector of x coordinates |
y |
vector of y coordinates |
r |
vactor of radii |
... |
additional graphic parameters will be passed through |
This function needs a previous plot where it adds the circles.
A. Gebhardt
x<-rnorm(10) y<-rnorm(10) r<-runif(10,0,0.5) plot(x,y, xlim=c(-3,3), ylim=c(-3,3), pch="+") circles(x,y,r)x<-rnorm(10) y<-rnorm(10) r<-runif(10,0,0.5) plot(x,y, xlim=c(-3,3), ylim=c(-3,3), pch="+") circles(x,y,r)
Sample data for the link{circumcircle} function.
circtest2 are points sampled from a circle with some jitter
added, i.e. they represent the most complicated case for the
link{circumcircle} function.
This function returns the circumcircle of a triangle and some additonal values used to determine them.
circum(x, y)circum(x, y)
x |
Vector of three elements, giving the x coordinatres of the triangle nodes. |
y |
Vector of three elements, giving the y coordinatres of the triangle nodes. |
This is an interface to the Fortran function CIRCUM found in TRIPACK.
x |
'x' coordinate of center |
y |
'y' coordinate of center |
radius |
circumcircle radius |
signed.area |
signed area of riangle (positive iff nodes are numbered counter clock wise) |
aspect.ratio |
ratio "radius of inscribed circle"/"radius of circumcircle", varies between 0 and 0.5 0 means collinear points, 0.5 equilateral trangle. |
This function is mainly intended to be used by circumcircle.
A. Gebhardt
https://math.fandom.com/wiki/Circumscribed_circle#Coordinates_of_circumcenter, visited march 2022.
circum(c(0,1,0),c(0,0,1)) tr <- list() tr$t1 <-list(x=c(0,1,0),y=c(0,0,1)) tr$t2 <-list(x=c(0.5,0.9,0.7),y=c(0.2,0.9,1)) tr$t3 <-list(x=c(0.05,0,0.3),y=c(0.2,0.7,0.1)) plot(0,0,type="n",xlim=c(-0.5,1.5),ylim=c(-0.5,1.5)) for(i in 1:3){ x <- tr[[i]]$x y <- tr[[i]]$y points(x,y,pch=c("1","2","3"),xlim=c(-0.5,1.5),ylim=c(-0.5,1.5)) cc =circum(x,y) lines(c(x,x[1]),c(y,y[1])) points(cc$x,cc$y) if(cc$signed.area<0) circles(cc$x,cc$y,cc$radius,col="blue",lty="dotted") else circles(cc$x,cc$y,cc$radius,col="red",lty="dotted") }circum(c(0,1,0),c(0,0,1)) tr <- list() tr$t1 <-list(x=c(0,1,0),y=c(0,0,1)) tr$t2 <-list(x=c(0.5,0.9,0.7),y=c(0.2,0.9,1)) tr$t3 <-list(x=c(0.05,0,0.3),y=c(0.2,0.7,0.1)) plot(0,0,type="n",xlim=c(-0.5,1.5),ylim=c(-0.5,1.5)) for(i in 1:3){ x <- tr[[i]]$x y <- tr[[i]]$y points(x,y,pch=c("1","2","3"),xlim=c(-0.5,1.5),ylim=c(-0.5,1.5)) cc =circum(x,y) lines(c(x,x[1]),c(y,y[1])) points(cc$x,cc$y) if(cc$signed.area<0) circles(cc$x,cc$y,cc$radius,col="blue",lty="dotted") else circles(cc$x,cc$y,cc$radius,col="red",lty="dotted") }
This function returns the (smallest) circumcircle of a set of n points
circumcircle(x, y = NULL, num.touch=2, plot = FALSE, debug = FALSE)circumcircle(x, y = NULL, num.touch=2, plot = FALSE, debug = FALSE)
x |
vector containing x coordinates of the data. If |
y |
vector containing y coordinates of the data. |
num.touch |
How often should the resulting circle touch the convex hull of the given points? default: 2 possible values: 2 or 3 Note: The circumcircle of a triangle is usually defined to touch
at 3 points, this function searches by default the minimum circle,
which may be only touching at 2 points. Set parameter
|
plot |
Logical, produce a simple plot of the result. default: |
debug |
Logical, more plots, only needed for debugging. default: |
This is a (naive implemented) algorithm which determines the smallest circumcircle of n points:
First step: Take the convex hull.
Second step: Determine two points on the convex hull with maximum distance for the diameter of the set.
Third step: Check if the circumcircle of these two points already contains all other points (of the convex hull and hence all other points).
If not or if 3 or more touching points are desired
(num.touch=3),
search a point with minimum enclosing circumcircle among the
remaining points of the convex hull.
If such a point cannot be found (e.g. for data(circtest2)),
search the remaining triangle combinations of points from the convex
hull until an enclosing circle with minimum radius is found.
The last search uses an upper and lower bound for the desired miniumum radius:
Any enclosing rectangle and its circumcircle gives an upper bound (the axis-parallel rectangle is used).
Half the diameter of the set from step 1 is a lower bound.
x |
'x' coordinate of circumcircle center |
y |
'y' coordinate of circumcircle center |
radius |
radius of circumcircle |
Albrecht Gebhardt
data(circtest) # smallest circle: circumcircle(circtest,num.touch=2,plot=TRUE) # smallest circle with maximum touching points (3): circumcircle(circtest,num.touch=3,plot=TRUE) # some stress test for this function, data(circtest2) # circtest2 was generated by: # 100 random points almost one a circle: # alpha <- runif(100,0,2*pi) # x <- cos(alpha) # y <- sin(alpha) # circtest2<-list(x=cos(alpha)+runif(100,0,0.1), # y=sin(alpha)+runif(100,0,0.1)) # circumcircle(circtest2,plot=TRUE)data(circtest) # smallest circle: circumcircle(circtest,num.touch=2,plot=TRUE) # smallest circle with maximum touching points (3): circumcircle(circtest,num.touch=3,plot=TRUE) # some stress test for this function, data(circtest2) # circtest2 was generated by: # 100 random points almost one a circle: # alpha <- runif(100,0,2*pi) # x <- cos(alpha) # y <- sin(alpha) # circtest2<-list(x=cos(alpha)+runif(100,0,0.1), # y=sin(alpha)+runif(100,0,0.1)) # circumcircle(circtest2,plot=TRUE)
Given a triangulation tri.obj of points in the plane, this
subroutine returns two vectors containing the coordinates
of the nodes on the boundary of the convex hull.
ConvexHull is an experimental C++ implementation of Grahams Scan
without previous triangulation, should be much faster.
convex.hull(tri.obj, plot.it=FALSE, add=FALSE,...) ConvexHull(x,y)convex.hull(tri.obj, plot.it=FALSE, add=FALSE,...) ConvexHull(x,y)
tri.obj |
object of class |
plot.it |
logical, if |
add |
logical. if |
... |
additional plot arguments |
x |
only for |
y |
only for |
x |
x coordinates of boundary nodes. |
y |
y coordinates of boundary nodes. |
In case that there are several collinear nodes on the convex hull
convex.hull will return them all while ConvexHull will only
give edge points.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht, print.triSht,
plot.triSht, summary.triSht, triangles.
## random points: rand.tr<-tri.mesh(runif(10),runif(10)) plot(rand.tr) rand.ch<-convex.hull(rand.tr, plot.it=TRUE, add=TRUE, col="red") ## use a part of the quakes data set: data(quakes) quakes.part<-quakes[(quakes[,1]<=-17 & quakes[,1]>=-19.0 & quakes[,2]<=182.0 & quakes[,2]>=180.0),] quakes.tri<-tri.mesh(quakes.part$lon, quakes.part$lat, duplicate="remove") plot(quakes.tri) convex.hull(quakes.tri, plot.it=TRUE, add=TRUE, col="red")## random points: rand.tr<-tri.mesh(runif(10),runif(10)) plot(rand.tr) rand.ch<-convex.hull(rand.tr, plot.it=TRUE, add=TRUE, col="red") ## use a part of the quakes data set: data(quakes) quakes.part<-quakes[(quakes[,1]<=-17 & quakes[,1]>=-19.0 & quakes[,2]<=182.0 & quakes[,2]>=180.0),] quakes.tri<-tri.mesh(quakes.part$lon, quakes.part$lat, duplicate="remove") plot(quakes.tri) convex.hull(quakes.tri, plot.it=TRUE, add=TRUE, col="red")
franke.data generates the test datasets from Franke, 1979, see references.
franke.data(fn = 1, ds = 1, data) franke.fn(x, y, fn = 1)franke.data(fn = 1, ds = 1, data) franke.fn(x, y, fn = 1)
fn |
function number, from 1 to 5. |
x |
'x' value |
y |
'y' value |
ds |
data set number, from 1 to 3. Dataset 1 consists of 100 points,
dataset 2 of 33 points and dataset 3 of 25 points scattered in the
square |
data |
A list of dataframes with 'x' and 'y' to choose from, dataset
|
These datasets are mentioned in Akima, (1996) as a testbed for the irregular scattered data interpolator.
Franke used the five functions:
and evaluated them on different more or less dense grids over .
A data frame with components
x |
'x' coordinate |
y |
'y' coordinate |
z |
'z' value |
The datasets have to be generated via franke.data before
use, the dataset franke only contains a list of 3 dataframes of
'x' and 'y' coordinates for the above mentioned irregular grids.
Do not forget to load the franke dataset first.
The 'x' and 'y' values have been taken from Akima (1996).
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
FRANKE, R., (1979). A critical comparison of some methods for interpolation of scattered data. Tech. Rep. NPS-53-79-003, Dept. of Mathematics, Naval Postgraduate School, Monterey, Calif.
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.
## generate Frankes data set for function 2 and dataset 3: data(franke) F23 <- franke.data(2,3,franke) str(F23)## generate Frankes data set for function 2 and dataset 3: data(franke) F23 <- franke.data(2,3,franke) str(F23)
Identify points in a plot of "x" with its
coordinates. The plot of "x" must be generated with plot.tri.
## S3 method for class 'triSht' identify(x,...)## S3 method for class 'triSht' identify(x,...)
x |
object of class |
... |
additional paramters for |
an integer vector containing the indexes of the identified points.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht, print.triSht, plot.triSht, summary.triSht
## Not run: data(franke) tr <- tri.mesh(franke$ds3$x, franke$ds3$y) plot(tr) identify(tr) ## End(Not run)## Not run: data(franke) tr <- tri.mesh(franke$ds3$x, franke$ds3$y) plot(tr) identify(tr) ## End(Not run)
This function implements bivariate interpolation for irregularly spaced input data. Piecewise linear (=barycentric interpolation), bilinear or bicubic spline interpolation according to Akimas method is applied.
interp(x, y = NULL, z, xo = seq(min(x), max(x), length = nx), yo = seq(min(y), max(y), length = ny), linear = (method == "linear"), extrap = FALSE, duplicate = "error", dupfun = NULL, nx = 40, ny = 40, input="points", output = "grid", method = "linear", deltri = "shull", h=0, kernel="gaussian", solver="QR", degree=3, baryweight=TRUE, autodegree=FALSE, adtol=0.1, smoothpde=FALSE, akimaweight=TRUE, nweight=25, na.rm=FALSE)interp(x, y = NULL, z, xo = seq(min(x), max(x), length = nx), yo = seq(min(y), max(y), length = ny), linear = (method == "linear"), extrap = FALSE, duplicate = "error", dupfun = NULL, nx = 40, ny = 40, input="points", output = "grid", method = "linear", deltri = "shull", h=0, kernel="gaussian", solver="QR", degree=3, baryweight=TRUE, autodegree=FALSE, adtol=0.1, smoothpde=FALSE, akimaweight=TRUE, nweight=25, na.rm=FALSE)
x |
vector of |
y |
vector of If left as NULL indicates that |
z |
vector of Missing values are not accepted by default, see parameter
|
xo |
If If |
yo |
If If |
input |
text, possible values are This is used to distinguish between regular and irregular gridded input data. |
output |
text, possible values are If In the case of |
linear |
logical, only for backward compatibility with Please use the new |
method |
text, possible methods are
|
extrap |
logical, indicates if extrapolation outside the convex hull is intended, this will not work for piecewise linear interpolation! |
duplicate |
character string indicating how to handle duplicate data points. Possible values are
|
dupfun |
a function, applied to duplicate points if
|
nx |
dimension of output grid in x direction |
ny |
dimension of output grid in y direction |
deltri |
triangulation method used, this argument may later be moved
into a control set together with others related to the spline
interpolation! Possible values are |
h |
bandwidth for partial derivatives estimation, compare |
kernel |
kernel for partial derivatives estimation, compare |
solver |
solver used in partial derivatives estimation, compare |
degree |
degree of local polynomial used for partial derivatives
estimation, compare |
baryweight |
calculate three partial derivatives estimators and return a barycentric weighted average. This increases the accuracy of Akima splines but the runtime is multplied by 3! |
autodegree |
try to reduce |
adtol |
tolerance used for autodegree |
smoothpde |
Use an averaged version of partial derivatives
estimates, by default simple average of Currently disabled by default (FALSE), underlying code still a bit experimental. |
akimaweight |
apply Akima weighting scheme on partial derivatives estimations instead of simply averaging |
nweight |
size of search neighbourhood for weighting scheme, default: 25 |
na.rm |
remove points where z= |
a list with 3 components:
x, y
|
If If |
z |
If If If the input was a |
Please note that this function tries to be a replacement for the interp() function from the akima package. So it should be call compatible for most applications. It also offers additional tuning parameters, usually the default settings will fit. Please be aware that these additional parameters may change in the future as they are still under development.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
Moebius, A. F. (1827) Der barymetrische Calcul. Verlag v. Johann Ambrosius Barth, Leipzig, https://books.google.at/books?id=eFPluv_UqFEC&hl=de&pg=PR1#v=onepage&q&f=false
Franke, R., (1979). A critical comparison of some methods for interpolation of scattered data. Tech. Rep. NPS-53-79-003, Dept. of Mathematics, Naval Postgraduate School, Monterey, Calif.
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.
### Use all datasets from Franke, 1979: data(franke) ## x-y irregular grid points: oldseed <- set.seed(42) ni <- 64 xi <- runif(ni,0,1) yi <- runif(ni,0,1) xyi <- cbind(xi,yi) ## linear interpolation fi <- franke.fn(xi,yi,1) IL <- interp(xi,yi,fi,nx=80,ny=80,method="linear") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IL$z,na.rm=TRUE),max(IL$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IL,breaks=breaks,col=colors,main="Franke function 1", sub=paste("linear interpolation, ", ni,"points")) contour(IL,add=TRUE,levels=breaks) points(xi,yi) ## spline interpolation fi <- franke.fn(xi,yi,1) IS <- interp(xi,yi,fi,method="akima", kernel="gaussian",solver="QR") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IS$z,na.rm=TRUE),max(IS$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IS,breaks=breaks,col=colors,main="Franke function 1", sub=paste("spline interpolation, ", ni,"points")) contour(IS,add=TRUE,levels=breaks) points(xi,yi) ## regular grid: nx <- 8; ny <- 8 xg<-seq(0,1,length=nx) yg<-seq(0,1,length=ny) xx <- t(matrix(rep(xg,ny),nx,ny)) yy <- matrix(rep(yg,nx),ny,nx) xyg<-expand.grid(xg,yg) ## linear interpolation fg <- outer(xg,yg,function(x,y)franke.fn(x,y,1)) IL <- interp(xg,yg,fg,input="grid",method="linear") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IL$z,na.rm=TRUE),max(IL$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IL,breaks=breaks,col=colors,main="Franke function 1", sub=paste("linear interpolation, ", nx,"x",ny,"points")) contour(IL,add=TRUE,levels=breaks) points(xx,yy) ## spline interpolation fg <- outer(xg,yg,function(x,y)franke.fn(x,y,1)) IS <- interp(xg,yg,fg,input="grid",method="akima", kernel="gaussian",solver="QR") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IS$z,na.rm=TRUE),max(IS$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IS,breaks=breaks,col=colors,main="Franke function 1", sub=paste("spline interpolation, ", nx,"x",ny,"points")) contour(IS,add=TRUE,levels=breaks) points(xx,yy) ## apply interp to sp data: require(sp) ## convert Akima data set to a sp object data(akima) asp <- SpatialPointsDataFrame(list(x=akima$x,y=akima$y), data = data.frame(z=akima$z)) spplot(asp,"z") ## linear interpolation spli <- interp(asp, z="z", method="linear") ## the result is again a SpatialPointsDataFrame: spplot(spli,"z") ## now with spline interpolation, slightly higher resolution spsi <- interp(asp, z="z", method="akima", nx=120, ny=120) spplot(spsi,"z") ## now sp grids: reuse stuff from above spgr <- SpatialPixelsDataFrame(list(x=c(xx),y=c(yy)), data=data.frame(z=c(fg))) spplot(spgr) ## linear interpolation spli <- interp(spgr, z="z", method="linear", input="grid") ## the result is again a SpatialPointsDataFrame: spplot(spli,"z") ## now with spline interpolation, slightly higher resolution spsi <- interp(spgr, z="z", method="akima", nx=240, ny=240) spplot(spsi,"z") set.seed(oldseed)### Use all datasets from Franke, 1979: data(franke) ## x-y irregular grid points: oldseed <- set.seed(42) ni <- 64 xi <- runif(ni,0,1) yi <- runif(ni,0,1) xyi <- cbind(xi,yi) ## linear interpolation fi <- franke.fn(xi,yi,1) IL <- interp(xi,yi,fi,nx=80,ny=80,method="linear") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IL$z,na.rm=TRUE),max(IL$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IL,breaks=breaks,col=colors,main="Franke function 1", sub=paste("linear interpolation, ", ni,"points")) contour(IL,add=TRUE,levels=breaks) points(xi,yi) ## spline interpolation fi <- franke.fn(xi,yi,1) IS <- interp(xi,yi,fi,method="akima", kernel="gaussian",solver="QR") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IS$z,na.rm=TRUE),max(IS$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IS,breaks=breaks,col=colors,main="Franke function 1", sub=paste("spline interpolation, ", ni,"points")) contour(IS,add=TRUE,levels=breaks) points(xi,yi) ## regular grid: nx <- 8; ny <- 8 xg<-seq(0,1,length=nx) yg<-seq(0,1,length=ny) xx <- t(matrix(rep(xg,ny),nx,ny)) yy <- matrix(rep(yg,nx),ny,nx) xyg<-expand.grid(xg,yg) ## linear interpolation fg <- outer(xg,yg,function(x,y)franke.fn(x,y,1)) IL <- interp(xg,yg,fg,input="grid",method="linear") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IL$z,na.rm=TRUE),max(IL$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IL,breaks=breaks,col=colors,main="Franke function 1", sub=paste("linear interpolation, ", nx,"x",ny,"points")) contour(IL,add=TRUE,levels=breaks) points(xx,yy) ## spline interpolation fg <- outer(xg,yg,function(x,y)franke.fn(x,y,1)) IS <- interp(xg,yg,fg,input="grid",method="akima", kernel="gaussian",solver="QR") ## prepare breaks and colors that match for image and contour: breaks <- pretty(seq(min(IS$z,na.rm=TRUE),max(IS$z,na.rm=TRUE),length=11)) db <- breaks[2]-breaks[1] nb <- length(breaks) breaks <- c(breaks[1]-db,breaks,breaks[nb]+db) colors <- terrain.colors(length(breaks)-1) image(IS,breaks=breaks,col=colors,main="Franke function 1", sub=paste("spline interpolation, ", nx,"x",ny,"points")) contour(IS,add=TRUE,levels=breaks) points(xx,yy) ## apply interp to sp data: require(sp) ## convert Akima data set to a sp object data(akima) asp <- SpatialPointsDataFrame(list(x=akima$x,y=akima$y), data = data.frame(z=akima$z)) spplot(asp,"z") ## linear interpolation spli <- interp(asp, z="z", method="linear") ## the result is again a SpatialPointsDataFrame: spplot(spli,"z") ## now with spline interpolation, slightly higher resolution spsi <- interp(asp, z="z", method="akima", nx=120, ny=120) spplot(spsi,"z") ## now sp grids: reuse stuff from above spgr <- SpatialPixelsDataFrame(list(x=c(xx),y=c(yy)), data=data.frame(z=c(fg))) spplot(spgr) ## linear interpolation spli <- interp(spgr, z="z", method="linear", input="grid") ## the result is again a SpatialPointsDataFrame: spplot(spli,"z") ## now with spline interpolation, slightly higher resolution spsi <- interp(spgr, z="z", method="akima", nx=240, ny=240) spplot(spsi,"z") set.seed(oldseed)
From an interp() result, produce a 3-column matrix
or data.frame cbind(x, y, z).
interp2xyz(al, data.frame = FALSE)interp2xyz(al, data.frame = FALSE)
al |
|
data.frame |
logical indicating if result should be data.frame or matrix (default). |
a matrix (or data.frame) with three columns, called
"x", "y", "z".
Martin Maechler, Jan.18, 2013
expand.grid() is the “essential ingredient” of
interp2xyz().
data(akima) ak.spl <- with(akima, interp(x, y, z, method = "akima")) str(ak.spl)# list (x[i], y[j], z = <matrix>[i,j]) ## Now transform to simple (x,y,z) matrix / data.frame : str(am <- interp2xyz(ak.spl)) str(ad <- interp2xyz(ak.spl, data.frame=TRUE)) ## and they are the same: stopifnot( am == ad | (is.na(am) & is.na(ad)) )data(akima) ak.spl <- with(akima, interp(x, y, z, method = "akima")) str(ak.spl)# list (x[i], y[j], z = <matrix>[i,j]) ## Now transform to simple (x,y,z) matrix / data.frame : str(am <- interp2xyz(ak.spl)) str(ad <- interp2xyz(ak.spl, data.frame=TRUE)) ## and they are the same: stopifnot( am == ad | (is.na(am) & is.na(ad)) )
This function implements bivariate interpolation onto a set of points for irregularly spaced input data.
This function is meant for backward compatibility to package
akima, please use interp with its output
argument set to "points" now. Especially newer options to the underlying
algorithm are only available there.
interpp(x, y = NULL, z, xo, yo = NULL, linear = TRUE, extrap = FALSE, duplicate = "error", dupfun = NULL, deltri = "shull")interpp(x, y = NULL, z, xo, yo = NULL, linear = TRUE, extrap = FALSE, duplicate = "error", dupfun = NULL, deltri = "shull")
x |
vector of x-coordinates of data points or a
|
y |
vector of y-coordinates of data points. Missing values are not accepted. If left as NULL indicates that |
z |
vector of z-coordinates of data points or a character variable
naming the variable of interest in the
Missing values are not accepted.
|
xo |
vector of x-coordinates of points at which to evaluate the interpolating
function. If |
yo |
vector of y-coordinates of points at which to evaluate the interpolating function. If operating on |
linear |
logical – indicating wether linear or spline interpolation should be used. |
extrap |
logical flag: should extrapolation be used outside of the convex hull determined by the data points? Not possible for linear interpolation. |
duplicate |
indicates how to handle duplicate data points. Possible values are
|
dupfun |
this function is applied to duplicate points if |
deltri |
triangulation method used, this argument will later be moved into a control set together with others related to the spline interpolation! |
a list with 3 components:
x, y
|
If If |
z |
If If If the input was a |
This is only a call wrapper meant for backward compatibility, see
interp for more details!
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
Moebius, A. F. (1827) Der barymetrische Calcul. Verlag v. Johann Ambrosius Barth, Leipzig, https://books.google.at/books?id=eFPluv_UqFEC&hl=de&pg=PR1#v=onepage&q&f=false
Franke, R., (1979). A critical comparison of some methods for interpolation of scattered data. Tech. Rep. NPS-53-79-003, Dept. of Mathematics, Naval Postgraduate School, Monterey, Calif.
### Use all datasets from Franke, 1979: ### calculate z at shifted original locations. data(franke) for(i in 1:5) for(j in 1:3){ FR <- franke.data(i,j,franke) IL <- with(FR, interpp(x,y,z,x+0.1,y+0.1,linear=TRUE)) str(IL) }### Use all datasets from Franke, 1979: ### calculate z at shifted original locations. data(franke) for(i in 1:5) for(j in 1:3){ FR <- franke.data(i,j,franke) IL <- with(FR, interpp(x,y,z,x+0.1,y+0.1,linear=TRUE)) str(IL) }
This function performs a local polynomial fit of up to order 3 to bivariate data. It returns estimated values of the regression function as well as estimated partial derivatives up to order 3. This access to the partial derivatives was the main intent for writing this code as there already many other local polynomial regression implementations in R.
locpoly(x, y, z, xo = seq(min(x), max(x), length = nx), yo = seq(min(y), max(y), length = ny), nx = 40, ny = 40, input = "points", output = "grid", h = 0, kernel = "gaussian", solver = "QR", degree = 3, pd = "")locpoly(x, y, z, xo = seq(min(x), max(x), length = nx), yo = seq(min(y), max(y), length = ny), nx = 40, ny = 40, input = "points", output = "grid", h = 0, kernel = "gaussian", solver = "QR", degree = 3, pd = "")
x |
vector of Missing values are not accepted. |
y |
vector of Missing values are not accepted. |
z |
vector of Missing values are not accepted.
|
xo |
If If |
yo |
If If |
input |
text, possible values are This is used to distinguish between regular and irregular gridded data. |
output |
text, possible values are If In the case of |
nx |
dimension of output grid in x direction |
ny |
dimension of output grid in y direction |
h |
bandwidth parameter, between 0 and 1. If a scalar is given it is interpreted as ratio applied to the dataset size to determine a local search neighbourhood, if set to 0 a minimum useful search neighbourhood is choosen (e.g. 10 points for a cubic trend function to determine all 10 parameters). If a vector of length 2 is given both components are interpreted as
ratio of the |
kernel |
Text value, implemented kernels are |
solver |
Text value, determines used solver in fastLM algorithm used by this code Possible values are |
degree |
Integer value, degree of polynomial trend, maximum allowed value is 3. |
pd |
Text value, determines which partial derivative should be returned,
possible values are |
If pd="all":
x |
|
y |
|
z |
estimates of |
zx |
estimates of |
zy |
estimates of |
zxx |
estimates of |
zxy |
estimates of |
zyy |
estimates of |
zxxx |
estimates of |
zxxy |
estimates of |
zxyy |
estimates of |
zyyy |
estimates of |
If pd!="all" only the elements x, y and the desired
derivative will be returned, e.g. zxy for pd="xy".
Function locpoly of package
KernSmooth performs a similar task for univariate data.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
Douglas Bates, Dirk Eddelbuettel (2013). Fast and Elegant Numerical Linear Algebra Using the RcppEigen Package. Journal of Statistical Software, 52(5), 1-24. URL http://www.jstatsoft.org/v52/i05/.
## choose a kernel knl <- "gaussian" ## choose global and local bandwidth bwg <- 0.25 # *100% means: percentage of x- y-range used bwl <- 0.1 # *100% means: percentage of data set (nearest neighbours) used ## a bivariate polynomial of degree 5: f <- function(x,y) 0.1+ 0.2*x-0.3*y+0.1*x*y+0.3*x^2*y-0.5*y^2*x+y^3*x^2+0.1*y^5 ## degree of model dg=3 ## part 1: ## regular gridded data: ng<- 11 # x/y size of a square data grid ## build and fill the grid with the theoretical values: xg<-seq(0,1,length=ng) yg<-seq(0,1,length=ng) # xg and yg as matrix matching fg nx <- length(xg) ny <- length(yg) xx <- t(matrix(rep(xg,ny),nx,ny)) yy <- matrix(rep(yg,nx),ny,nx) fg <- outer(xg,yg,f) ## local polynomial estimate ## global bw: ttg <- system.time(pdg <- locpoly(xg,yg,fg, input="grid", pd="all", h=c(bwg,bwg), solver="QR", degree=dg, kernel=knl)) ## time used: ttg ## local bw: ttl <- system.time(pdl <- locpoly(xg,yg,fg, input="grid", pd="all", h=bwl, solver="QR", degree=dg, kernel=knl)) ## time used: ttl image(pdl$x,pdl$y,pdl$z,main="f and its estimated first partial derivatives", sub="colors: f, dotted: df/dx, dashed: df/dy") contour(pdl$x,pdl$y,pdl$zx,add=TRUE,lty="dotted") contour(pdl$x,pdl$y,pdl$zy,add=TRUE,lty="dashed") points(xx,yy,pch=".") ## part 2: ## irregular data, ## results will not be as good as with the regular 21*21=231 points. nd<- 121 # size of data set ## random irregular data oldseed <- set.seed(42) x<-runif(ng) y<-runif(ng) set.seed(oldseed) z <- f(x,y) ## global bw: ttg <- system.time(pdg <- interp::locpoly(x,y,z, xg,yg, pd="all", h=c(bwg,bwg), solver="QR", degree=dg,kernel=knl)) ttg ## local bw: ttl <- system.time(pdl <- interp::locpoly(x,y,z, xg,yg, pd="all", h=bwl, solver="QR", degree=dg,kernel=knl)) ttl image(pdl$x,pdl$y,pdl$z,main="f and its estimated first partial derivatives", sub="colors: f, dotted: df/dx, dashed: df/dy") contour(pdl$x,pdl$y,pdl$zx,add=TRUE,lty="dotted") contour(pdl$x,pdl$y,pdl$zy,add=TRUE,lty="dashed") points(x,y,pch=".")## choose a kernel knl <- "gaussian" ## choose global and local bandwidth bwg <- 0.25 # *100% means: percentage of x- y-range used bwl <- 0.1 # *100% means: percentage of data set (nearest neighbours) used ## a bivariate polynomial of degree 5: f <- function(x,y) 0.1+ 0.2*x-0.3*y+0.1*x*y+0.3*x^2*y-0.5*y^2*x+y^3*x^2+0.1*y^5 ## degree of model dg=3 ## part 1: ## regular gridded data: ng<- 11 # x/y size of a square data grid ## build and fill the grid with the theoretical values: xg<-seq(0,1,length=ng) yg<-seq(0,1,length=ng) # xg and yg as matrix matching fg nx <- length(xg) ny <- length(yg) xx <- t(matrix(rep(xg,ny),nx,ny)) yy <- matrix(rep(yg,nx),ny,nx) fg <- outer(xg,yg,f) ## local polynomial estimate ## global bw: ttg <- system.time(pdg <- locpoly(xg,yg,fg, input="grid", pd="all", h=c(bwg,bwg), solver="QR", degree=dg, kernel=knl)) ## time used: ttg ## local bw: ttl <- system.time(pdl <- locpoly(xg,yg,fg, input="grid", pd="all", h=bwl, solver="QR", degree=dg, kernel=knl)) ## time used: ttl image(pdl$x,pdl$y,pdl$z,main="f and its estimated first partial derivatives", sub="colors: f, dotted: df/dx, dashed: df/dy") contour(pdl$x,pdl$y,pdl$zx,add=TRUE,lty="dotted") contour(pdl$x,pdl$y,pdl$zy,add=TRUE,lty="dashed") points(xx,yy,pch=".") ## part 2: ## irregular data, ## results will not be as good as with the regular 21*21=231 points. nd<- 121 # size of data set ## random irregular data oldseed <- set.seed(42) x<-runif(ng) y<-runif(ng) set.seed(oldseed) z <- f(x,y) ## global bw: ttg <- system.time(pdg <- interp::locpoly(x,y,z, xg,yg, pd="all", h=c(bwg,bwg), solver="QR", degree=dg,kernel=knl)) ttg ## local bw: ttl <- system.time(pdl <- interp::locpoly(x,y,z, xg,yg, pd="all", h=bwl, solver="QR", degree=dg,kernel=knl)) ttl image(pdl$x,pdl$y,pdl$z,main="f and its estimated first partial derivatives", sub="colors: f, dotted: df/dx, dashed: df/dy") contour(pdl$x,pdl$y,pdl$zx,add=TRUE,lty="dotted") contour(pdl$x,pdl$y,pdl$zy,add=TRUE,lty="dashed") points(x,y,pch=".")
This function can be used to generate nearest neighbour information for a set of 2D data points.
nearest.neighbours(x, y)nearest.neighbours(x, y)
x |
vector containing |
y |
vector containing |
The C++ implementation of this function is used inside the
locpoly and interp functions.
A list with two components
index |
A matrix with one row per data point. Each row contains the indices
of the nearest neigbours to the point associated with this row,
currently the point itself is also listed in the first row, so
this matrix is of dimension |
dist |
A matrix containing the distances according to the neigbours listed
in component |
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
data(franke) ## use only a small subset fd <- franke$ds1[1:5,] nearest.neighbours(fd$x,fd$y)data(franke) ## use only a small subset fd <- franke$ds1[1:5,] nearest.neighbours(fd$x,fd$y)
Extract a list of neighbours from a triangulation or voronoi object
neighbours(obj)neighbours(obj)
obj |
object of class |
nested list of neighbours per point
A. Gebhardt
triSht, print.triSht, plot.triSht, summary.triSht, triangles
data(tritest) tritest.tr<-tri.mesh(tritest$x,tritest$y) tritest.nb<-neighbours(tritest.tr)data(tritest) tritest.tr<-tri.mesh(tritest$x,tritest$y) tritest.nb<-neighbours(tritest.tr)
A simple test function to determine the position of one (or more) points relative to a vector spanned by two points.
on(x1, y1, x2, y2, x0, y0, eps = 1e-16) left(x1, y1, x2, y2, x0, y0, eps = 1e-16)on(x1, y1, x2, y2, x0, y0, eps = 1e-16) left(x1, y1, x2, y2, x0, y0, eps = 1e-16)
x1 |
|
y1 |
|
x2 |
|
y2 |
|
x0 |
vector of |
y0 |
vector of |
eps |
tolerance for checking if |
logical vector with the results of the test.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
in.convex.hull, on.convex.hull.
y <- x <- c(0,1) ## should be TRUE on(x[1],y[1],x[2],y[2],0.5,0.5) ## note the default setting of eps leading to on(x[1],y[1],x[2],y[2],0.5,0.50000000000000001) ## also be TRUE ## should be TRUE left(x[1],y[1],x[2],y[2],0.5,0.6) ## note the default setting of eps leading to left(x[1],y[1],x[2],y[2],0.5,0.50000000000000001) ## already resulting to FALSEy <- x <- c(0,1) ## should be TRUE on(x[1],y[1],x[2],y[2],0.5,0.5) ## note the default setting of eps leading to on(x[1],y[1],x[2],y[2],0.5,0.50000000000000001) ## also be TRUE ## should be TRUE left(x[1],y[1],x[2],y[2],0.5,0.6) ## note the default setting of eps leading to left(x[1],y[1],x[2],y[2],0.5,0.50000000000000001) ## already resulting to FALSE
Given a triangulation object tri.obj of points in the plane, this
subroutine returns a logical vector indicating if the points
lay on or in the convex hull of tri.obj.
on.convex.hull(tri.obj, x, y, eps=1E-16) in.convex.hull(tri.obj, x, y, eps=1E-16, strict=TRUE)on.convex.hull(tri.obj, x, y, eps=1E-16) in.convex.hull(tri.obj, x, y, eps=1E-16, strict=TRUE)
tri.obj |
object of class |
x |
vector of |
y |
vector of |
eps |
accuracy for checking the condition |
strict |
logical, default |
Logical vector.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht, print.triSht, plot.triSht,
summary.triSht, triangles,
convex.hull.
# use a part of the quakes data set: data(quakes) quakes.part<-quakes[(quakes[,1]<=-10.78 & quakes[,1]>=-19.4 & quakes[,2]<=182.29 & quakes[,2]>=165.77),] q.tri<-tri.mesh(quakes.part$lon, quakes.part$lat, duplicate="remove") on.convex.hull(q.tri,quakes.part$lon[1:20],quakes.part$lat[1:20]) # Check with part of data set: # Note that points on the hull (see above) get marked FALSE below: in.convex.hull(q.tri,quakes.part$lon[1:20],quakes.part$lat[1:20]) # If points both on the hull and in the interior of the hull are meant # disable strict mode: in.convex.hull(q.tri,quakes.part$lon[1:20],quakes.part$lat[1:20],strict=FALSE) # something completely outside: in.convex.hull(q.tri,c(170,180),c(-20,-10))# use a part of the quakes data set: data(quakes) quakes.part<-quakes[(quakes[,1]<=-10.78 & quakes[,1]>=-19.4 & quakes[,2]<=182.29 & quakes[,2]>=165.77),] q.tri<-tri.mesh(quakes.part$lon, quakes.part$lat, duplicate="remove") on.convex.hull(q.tri,quakes.part$lon[1:20],quakes.part$lat[1:20]) # Check with part of data set: # Note that points on the hull (see above) get marked FALSE below: in.convex.hull(q.tri,quakes.part$lon[1:20],quakes.part$lat[1:20]) # If points both on the hull and in the interior of the hull are meant # disable strict mode: in.convex.hull(q.tri,quakes.part$lon[1:20],quakes.part$lat[1:20],strict=FALSE) # something completely outside: in.convex.hull(q.tri,c(170,180),c(-20,-10))
This version of outer evaluates FUN
only on that part of the grid times that is enclosed within
the convex hull of the points .
This can be useful for spatial estimation if no extrapolation is wanted.
outer.convhull(cx,cy,px,py,FUN,duplicate="remove",...)outer.convhull(cx,cy,px,py,FUN,duplicate="remove",...)
cx |
x cordinates of grid |
cy |
y cordinates of grid |
px |
vector of x coordinates of points |
py |
vector of y coordinates of points |
FUN |
function to be evaluated over the grid |
duplicate |
indicates what to do with duplicate
|
... |
additional arguments for |
Matrix with values of FUN (NAs if outside the
convex hull).
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
x<-runif(20) y<-runif(20) z<-runif(20) z.lm<-lm(z~x+y) f.pred<-function(x,y) {predict(z.lm,data.frame(x=as.vector(x),y=as.vector(y)))} xg<-seq(0,1,0.05) yg<-seq(0,1,0.05) image(xg,yg,outer.convhull(xg,yg,x,y,f.pred)) points(x,y)x<-runif(20) y<-runif(20) z<-runif(20) z.lm<-lm(z~x+y) f.pred<-function(x,y) {predict(z.lm,data.frame(x=as.vector(x),y=as.vector(y)))} xg<-seq(0,1,0.05) yg<-seq(0,1,0.05) image(xg,yg,outer.convhull(xg,yg,x,y,f.pred)) points(x,y)
plots the triangulation object "x"
## S3 method for class 'triSht' plot(x, add = FALSE, xlim = range(x$x), ylim = range(x$y), do.points = TRUE, do.labels = FALSE, isometric = TRUE, do.circumcircles = FALSE, segment.lty = "dashed", circle.lty = "dotted", ...)## S3 method for class 'triSht' plot(x, add = FALSE, xlim = range(x$x), ylim = range(x$y), do.points = TRUE, do.labels = FALSE, isometric = TRUE, do.circumcircles = FALSE, segment.lty = "dashed", circle.lty = "dotted", ...)
x |
object of class |
add |
logical, if |
do.points |
logical, indicates if points should be
plotted. (default |
do.labels |
logical, indicates if points should be labelled.
(default |
xlim, ylim
|
x/y ranges for plot |
isometric |
generate an isometric plot (default |
do.circumcircles |
logical, indicates if circumcircles should be
plotted (default |
segment.lty |
line type for triangulation segments |
circle.lty |
line type for circumcircles |
... |
additional plot parameters |
None
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht, print.triSht,
summary.triSht
## random points plot(tri.mesh(rpois(100,lambda=20),rpois(100,lambda=20),duplicate="remove")) ## use a part of the quakes data set: data(quakes) quakes.part<-quakes[(quakes[,1]<=-10.78 & quakes[,1]>=-19.4 & quakes[,2]<=182.29 & quakes[,2]>=165.77),] quakes.tri<-tri.mesh(quakes.part$lon, quakes.part$lat, duplicate="remove") plot(quakes.tri) ## use the whole quakes data set ## (will not work with standard memory settings, hence commented out) ## plot(tri.mesh(quakes$lon, quakes$lat, duplicate="remove"), do.points=F)## random points plot(tri.mesh(rpois(100,lambda=20),rpois(100,lambda=20),duplicate="remove")) ## use a part of the quakes data set: data(quakes) quakes.part<-quakes[(quakes[,1]<=-10.78 & quakes[,1]>=-19.4 & quakes[,2]<=182.29 & quakes[,2]>=165.77),] quakes.tri<-tri.mesh(quakes.part$lon, quakes.part$lat, duplicate="remove") plot(quakes.tri) ## use the whole quakes data set ## (will not work with standard memory settings, hence commented out) ## plot(tri.mesh(quakes$lon, quakes$lat, duplicate="remove"), do.points=F)
Plots the mosaic "x".
Dashed lines are used for outer tiles of the mosaic.
## S3 method for class 'voronoi' plot(x,add=FALSE, xlim=c(min(x$tri$x)- 0.1*diff(range(x$tri$x)), max(x$tri$x)+ 0.1*diff(range(x$tri$x))), ylim=c(min(x$tri$y)- 0.1*diff(range(x$tri$y)), max(x$tri$y)+ 0.1*diff(range(x$tri$y))), all=FALSE, do.points=TRUE, main="Voronoi mosaic", sub=deparse(substitute(x)), isometric=TRUE, ...)## S3 method for class 'voronoi' plot(x,add=FALSE, xlim=c(min(x$tri$x)- 0.1*diff(range(x$tri$x)), max(x$tri$x)+ 0.1*diff(range(x$tri$x))), ylim=c(min(x$tri$y)- 0.1*diff(range(x$tri$y)), max(x$tri$y)+ 0.1*diff(range(x$tri$y))), all=FALSE, do.points=TRUE, main="Voronoi mosaic", sub=deparse(substitute(x)), isometric=TRUE, ...)
x |
object of class |
add |
logical, if |
xlim |
x plot ranges, by default modified to hide dummy points outside of the plot |
ylim |
y plot ranges, by default modified to hide dummy points outside of the plot |
all |
show all (including dummy points in the plot |
do.points |
logical, indicates if points should be plotted. |
main |
plot title |
sub |
plot subtitle |
isometric |
generate an isometric plot (default |
... |
additional plot parameters |
None
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
voronoi, print.voronoi,
summary.voronoi, plot.voronoi.polygons
data(franke) tr <- tri.mesh(franke$ds3) vr <- voronoi.mosaic(tr) plot(tr) plot(vr,add=TRUE)data(franke) tr <- tri.mesh(franke$ds3) vr <- voronoi.mosaic(tr) plot(tr) plot(vr,add=TRUE)
plots an voronoi.polygons object
## S3 method for class 'voronoi.polygons' plot(x, which, color=TRUE, isometric=TRUE, ...)## S3 method for class 'voronoi.polygons' plot(x, which, color=TRUE, isometric=TRUE, ...)
x |
object of class |
which |
index vector selecting which polygons to plot |
color |
logical, determines if plot should be colored, default:
|
isometric |
generate an isometric plot (default |
... |
additional plot arguments |
A. Gebhardt
data(franke) fd3 <- franke$ds3 fd3.vm <- voronoi.mosaic(fd3$x,fd3$y) fd3.vp <- voronoi.polygons(fd3.vm) plot(fd3.vp) plot(fd3.vp,which=c(3,4,6,10))data(franke) fd3 <- franke$ds3 fd3.vm <- voronoi.mosaic(fd3$x,fd3$y) fd3.vp <- voronoi.polygons(fd3.vm) plot(fd3.vp) plot(fd3.vp,which=c(3,4,6,10))
Prints some information about tri.obj
## S3 method for class 'summary.triSht' print(x, ...)## S3 method for class 'summary.triSht' print(x, ...)
x |
object of class |
... |
additional paramters for |
None
This function is meant as replacement for the function of same name in
package tripack.
The only difference is that no constraints are possible with
triSht objects of package interp.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht,tri.mesh,
print.triSht, plot.triSht,
summary.triSht.
Prints some information about object x
## S3 method for class 'summary.voronoi' print(x, ...)## S3 method for class 'summary.voronoi' print(x, ...)
x |
object of class |
... |
additional paramters for |
None
This function is meant as replacement for the function of same name in
package tripack and should be fully backward compatible.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
voronoi,voronoi.mosaic,
print.voronoi, plot.voronoi,
summary.voronoi.
prints a adjacency list of "x"
## S3 method for class 'triSht' print(x,...)## S3 method for class 'triSht' print(x,...)
x |
object of class |
... |
additional paramters for |
None
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht,
plot.triSht,
summary.triSht
prints a summary of "x"
## S3 method for class 'voronoi' print(x,...)## S3 method for class 'voronoi' print(x,...)
x |
object of class |
... |
additional paramters for |
None
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
voronoi,
plot.voronoi,
summary.voronoi
Returns some information (number of nodes, triangles, arcs)
about object.
## S3 method for class 'triSht' summary(object,...)## S3 method for class 'triSht' summary(object,...)
object |
object of class |
... |
additional paramters for |
An object of class "summary.triSht", to be printed by
print.summary.triSht.
It contains the number of nodes (n), of arcs (na), of
boundary nodes (nb) and triangles (nt).
This function is meant as replacement for the function of same name in
package tripack.
The only difference is that no constraints are possible with
triSht objects of package interp.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht, print.triSht, plot.triSht,
print.summary.triSht.
Returns some information about object
## S3 method for class 'voronoi' summary(object,...)## S3 method for class 'voronoi' summary(object,...)
object |
object of class |
... |
additional parameters for |
Object of class "summary.voronoi".
It contains the number of nodes (nn) and dummy nodes (nd).
This function is meant as replacement for the function of same name in
package tripack and should be fully backward compatible.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
voronoi,voronoi.mosaic,
print.voronoi, plot.voronoi,
print.summary.voronoi.
This subroutine locates a point relative to a
triangulation created by tri.mesh. If is
contained in a triangle, the three vertex indexes are
returned. Otherwise, the indexes of the rightmost and
leftmost visible boundary nodes are returned.
tri.find(tri.obj,x,y)tri.find(tri.obj,x,y)
tri.obj |
an triangulation object of class |
x |
x-coordinate of the point |
y |
y-coordinate of the point |
A list with elements i1,i2,i3 containing nodal
indexes, in counterclockwise order, of the vertices of a triangle
containing . tr contains the triangle index and
bc contains the barycentric coordinates
of w.r.t. the found triangle.
If is not contained in the
convex hull of the nodes this indices are 0 (bc is meaningless then).
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht, print.triSht, plot.triSht,
summary.triSht, triangles,
convex.hull
data(franke) tr<-tri.mesh(franke$ds3$x,franke$ds3$y) plot(tr) pnt<-list(x=0.3,y=0.4) triangle.with.pnt<-tri.find(tr,pnt$x,pnt$y) attach(triangle.with.pnt) lines(franke$ds3$x[c(i1,i2,i3,i1)],franke$ds3$y[c(i1,i2,i3,i1)],col="red") points(pnt$x,pnt$y)data(franke) tr<-tri.mesh(franke$ds3$x,franke$ds3$y) plot(tr) pnt<-list(x=0.3,y=0.4) triangle.with.pnt<-tri.find(tr,pnt$x,pnt$y) attach(triangle.with.pnt) lines(franke$ds3$x[c(i1,i2,i3,i1)],franke$ds3$y[c(i1,i2,i3,i1)],col="red") points(pnt$x,pnt$y)
This function generates a Delaunay triangulation of arbitrarily distributed points in the plane. The resulting object can be printed or plotted, some additional functions can extract details from it like the list of triangles, arcs or the convex hull.
tri.mesh(x, y = NULL, duplicate = "error", jitter = FALSE)tri.mesh(x, y = NULL, duplicate = "error", jitter = FALSE)
x |
vector containing |
y |
vector containing |
duplicate |
flag indicating how to handle duplicate elements. Possible values are:
|
jitter |
logical, adds some jitter to both coordinates as this can
help in situations with too much colinearity. Default is |
This function creates a Delaunay triangulation of a set of arbitrarily distributed points in the plane referred to as nodes.
The Delaunay triangulation is defined as a set of triangles with the following five properties:
The triangle vertices are nodes.
No triangle contains a node other than its vertices.
The interiors of the triangles are pairwise disjoint.
The union of triangles is the convex hull of the set of nodes (the smallest convex set which contains the nodes).
The interior of the circumcircle of each triangle contains no node.
The first four properties define a triangulation, and the last property results in a triangulation which is as close as possible to equiangular in a certain sense and which is uniquely defined unless four or more nodes lie on a common circle. This property makes the triangulation well-suited for solving closest point problems and for triangle-based interpolation.
This triangulation is based on the s-hull algorithm by David Sinclair. It consist of two steps:
Create an initial non-overlapping triangulation from the radially sorted nodes (w.r.t to an arbitrary first node). Starting from a first triangle built from the first node and its nearest neigbours this is done by adding triangles from the next node (in the sense of distance to the first node) to the hull of the actual triangulation visible from this node (sweep hull step).
Apply triange flipping to each pair of triangles sharing a border until condition 5 holds (Cline-Renka test).
This algorithm has complexicity .
an object of class "triSht", see triSht.
This function is meant as a replacement for function
tri.mesh from package tripack.
Please note that the underlying algorithm changed from Renka's method
to Sinclair's sweep hull method. Delaunay triangulations are unique if
no four or more points exist which share the same
circumcircle. Otherwise several solutions are available and different
algorithms will give different results. This especially holds for
regular grids, where in the case of rectangular gridded points each
grid cell can be triangulated in two different ways.
The arguments are backward compatible, but the returned object
is not compatible with package tripack (it
provides a tri object type)! But you
can apply methods with same names to the object returned in package
interp which is of type triSht, so you can reuse
your old code but you cannot reuse your old saved workspace.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
B. Delaunay, Sur la sphere vide. A la memoire de Georges Voronoi, Bulletin de l'Academie des Sciences de l'URSS. Classe des sciences mathematiques et na, 1934, no. 6, p. 793–800
D. A. Sinclair, S-Hull: A Fast Radial Sweep-Hull Routine for Delaunay Triangulation. https://arxiv.org/pdf/1604.01428.pdf, 2016.
triSht, print.triSht, plot.triSht,
summary.triSht, triangles,
convex.hull, arcs.
## use Frankes datasets: data(franke) tr1 <- tri.mesh(franke$ds3$x, franke$ds3$y) tr1 tr2 <- tri.mesh(franke$ds2) summary(tr2)## use Frankes datasets: data(franke) tr1 <- tri.mesh(franke$ds3$x, franke$ds3$y) tr1 tr2 <- tri.mesh(franke$ds2) summary(tr2)
This function extracts a list of triangles
from an triangulation object created by tri.mesh.
triangles(tri.obj)triangles(tri.obj)
tri.obj |
object of class |
The vertices in the returned matrix (let's denote it with
retval) are ordered counterclockwise. The columns tr
and arc, index the triangle and arc, respectively,
which are opposite (not shared by) node node, with
tri if arc indexes a boundary arc. Vertex indexes
range from 1 to , the number of nodes, triangle indexes from 0
to , and arc indexes from 1 to .
A matrix with columns node1, node2, node3,
representing the vertex nodal indexes,
tr1, tr2, tr3, representing neighboring triangle
indexes and arc1, arc2, arc3 reresenting arc
indexes.
Each row represents one triangle.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
triSht, print.triSht,
plot.triSht, summary.triSht,
triangles
# use the smallest Franke data set data(franke) fr3.tr<-tri.mesh(franke$ds3$x, franke$ds3$y) triangles(fr3.tr)# use the smallest Franke data set data(franke) fr3.tr<-tri.mesh(franke$ds3$x, franke$ds3$y) triangles(fr3.tr)
R object that represents the triangulation of a set of 2D points,
generated by tri.mesh.
n |
Number of nodes |
x |
|
y |
|
nt |
number of triangles |
trlist |
Matrix of indices which defines the triangulation, each row corresponds to a triangle. Columns Columns Columns |
cclist |
Matrix describing the circumcircles and triangles. Columns
The radius of the inscribed circle can be get via
|
nchull |
number of points on the convex hull |
chull |
A vector containing the indices of nodes forming the convec hull (in counterclockwise ordering). |
narcs |
number of arcs forming the triangulation |
arcs |
A matrix with node indices describing the arcs, contains
two columns |
call |
call, which generated this object |
This object is not backward compatible with tri objects generated
from package tripack but the functions and methods are! So you
have to regenerate these objects and then you can continue to use the
same calls as before.
The only difference is that no constraints to the triangulation are
possible in package interp.
Function triSht2tri provides an option to convert this object into
the older form from package tripack, but it will not generate exact
copies as if the object would have been created with tripack::tri.mesh!
The old data structure consists of three lists describing adjacency lists
of triangulation nodes in counterclockwise order, the translation function
only genrates such a valid (but not unique) description.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
tri.mesh, print.triSht,triSht2tri,
plot.triSht, summary.triSht
This function converts triSht objects (from this package) to tri objects (from tripack package).
triSht2tri(t.triSht)triSht2tri(t.triSht)
t.triSht |
a class |
A class tri object, see tripack package.
The converted objects are not fully compatible with tripack functions. Basic stuff (printing, plotting) works, tripack::triangles e.g. does not work.
Voronoi functions from package tripack are working correctly with translated objects.
A. Gebhardt
A very simply set set of points to test the tripack functions, taken
from the FORTRAN original. tritest2 is a slight modification by
adding runif(,-0.1,0.1) random numbers to the coordinates.
R. J. Renka (1996). Algorithm 751: TRIPACK: a constrained two-dimensional Delaunay triangulation package. ACM Transactions on Mathematical Software. 22, 1-8.
A voronoi object is created with voronoi.mosaic
x, y
|
x and y coordinates of nodes of the voronoi mosaic. Each node is a circumcircle center of some triangle from the Delaunay triangulation. |
node |
logical vector, indicating real nodes of the voronoi mosaic. These nodes are the centers of circumcircles of triangles with positive area of the delaunay triangulation. If |
n1, n2, n3
|
indices of neighbour nodes. Negative indices indicate dummy points as neighbours. |
tri |
triangulation object, see |
area |
area of triangle |
ratio |
aspect ratio (inscribed radius/circumradius) of triangle
|
radius |
circumradius of triangle i. |
dummy.x, dummy.y
|
x and y coordinates of dummy points. They are used for plotting of unbounded tiles. |
This version of voronoi object is generated from the
tri.mesh function from package interp. That's the only
difference to voronoi objects generated with package
tripack.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
Computes the area of each Voronoi polygon.
For some sites at the edge of the region, the Voronoi polygon is not
bounded, and so the area of those sites cannot be calculated, and hence
will be NA.
voronoi.area(voronoi.obj)voronoi.area(voronoi.obj)
voronoi.obj |
object of class |
A vector of polygon areas.
S. J. Eglen
voronoi.mosaic,voronoi.polygons,
data(franke) fd3 <- franke$ds3 fd3.vm <- voronoi.mosaic(fd3$x,fd3$y) fd3.vm.areas <- voronoi.area(fd3.vm) plot(fd3.vm) text(fd3$x, fd3$y, round(fd3.vm.areas,5))data(franke) fd3 <- franke$ds3 fd3.vm <- voronoi.mosaic(fd3$x,fd3$y) fd3.vm.areas <- voronoi.area(fd3.vm) plot(fd3.vm) text(fd3$x, fd3$y, round(fd3.vm.areas,5))
Find the sites in the Voronoi tesselation that lie at the edge of the region. A site is at the edge if any of the vertices of its Voronoi polygon lie outside the rectangle with corners (xmin,ymin) and (xmax,ymax).
voronoi.findrejectsites(voronoi.obj, xmin, xmax, ymin, ymax)voronoi.findrejectsites(voronoi.obj, xmin, xmax, ymin, ymax)
voronoi.obj |
object of class |
xmin |
minimum x-coordinate of sites in the region |
xmax |
maximum x-coordinate of sites in the region |
ymin |
minimum y-coordinate of sites in the region |
ymax |
maximum y-coordinate of sites in the region |
A logical vector of the same length as the number of sites. If the site is a reject, the corresponding element of the vector is set to TRUE.
S. J. Eglen
This function creates a Voronoi mosaic out of a given set of
arbitraryly located points in the plane. Each cell of a voronoi
mosaic is associated with a data point and contains all points
closest to this data point.
voronoi.mosaic(x, y = NULL, duplicate = "error")voronoi.mosaic(x, y = NULL, duplicate = "error")
x |
vector containing
|
y |
vector containing |
duplicate |
flag indicating how to handle duplicate elements. Possible values are:
|
The function creates first a Delaunay triangulation (if not already given), extracts the circumcircle centers of these triangles, and then connects these points according to the neighbourhood relations between the triangles.
An object of class voronoi.
This function is meant as a replacement for function
voronoi.mosaic from package tripack.
Please note that the underlying triangulation uses a
different algorithm, see tri.mesh. Contrary to
tri.mesh this should not affect the result for non
unique triangulations e.g. on regular grids as the voronoi mosaic in
this case will still be unique.
The arguments are backward compatible, even the returned object should be
compatible with functions from package tripack.
Albrecht Gebhardt <[email protected]>, Roger Bivand <[email protected]>
G. Voronoi, Nouvelles applications des parametres continus a la theorie des formes quadratiques. Deuxieme memoire. Recherches sur les parallelloedres primitifs, Journal fuer die reine und angewandte Mathematik, 1908, vol 134, p. 198-287
voronoi,voronoi.mosaic,
print.voronoi, plot.voronoi
data(franke) fd <- franke$ds3 vr <- voronoi.mosaic(fd$x, fd$y) summary(vr)data(franke) fd <- franke$ds3 vr <- voronoi.mosaic(fd$x, fd$y) summary(vr)
This functions extracts polygons from a voronoi.mosaic object.
voronoi.polygons(voronoi.obj)voronoi.polygons(voronoi.obj)
voronoi.obj |
object of class |
Returns an object of class voronoi.polygons with unamed list
elements for each polygon. These list
elements are matrices with columns x and y.
Unbounded polygons along the border are represented by NULL
instead of a matrix.
Denis White
plot.voronoi.polygons,voronoi.mosaic
data(franke) fd3 <- franke$ds3 fd3.vm <- voronoi.mosaic(fd3$x,fd3$y) fd3.vp <- voronoi.polygons(fd3.vm) fd3.vpdata(franke) fd3 <- franke$ds3 fd3.vm <- voronoi.mosaic(fd3$x,fd3$y) fd3.vp <- voronoi.polygons(fd3.vm) fd3.vp