Skip to content

GitLab

Commit b9de5f5f authored by agebhard's avatar agebhard
Browse files

* re-add Akimas old code temporarily 

* declare interp.new and interp.old deprecated
parent fdd15e4a
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ b9de5f5f
Package: akima
Version: 0.6-1
Version: 0.6-2
Date: 2016-12-16
Title: Interpolation of Irregularly and Regularly Spaced Data
Authors@R: c(person("Hiroshi", "Akima", role = c("aut", "cph"),
comment = "Fortran code (TOMS 526, 761, 697 and 433)"),
comment = "Fortran code (TOMS 760, 761, 697 and 433)"),
person("Albrecht", "Gebhardt", role = c("aut", "cre", "cph"),
email = "albrecht.gebhardt@aau.at",
comment = "R port (interp* functions), bicubic* functions"),
......@@ -12,7 +12,7 @@ Authors@R: c(person("Hiroshi", "Akima", role = c("aut", "cph"),
person("Martin", "Maechler", email="maechler@stat.math.ethz.ch",
role = c("ctb", "cph"), comment = "interp2xyz function + enhancements"),
person("YYYY Association for Computing Machinery, Inc.",
role = c("cph"), comment = "covers code from TOMS 526, 761, 697 and 433")
role = c("cph"), comment = "covers code from TOMS 760, 761, 697 and 433")
)
Maintainer: Albrecht Gebhardt <albrecht.gebhardt@aau.at>
Description: Several cubic spline interpolation methods of H. Akima for irregular and
......
NAMESPACE
View file @ b9de5f5f
useDynLib(akima)
export(aspline,interp,interpp,interp2xyz,bicubic,bicubic.grid,bilinear,bilinear.grid,franke.fn,franke.data)
export(aspline,interp,interpp,interp.old,interp.new,interpp.old,interpp.new,
interp2xyz,
bicubic,bicubic.grid,bilinear,bilinear.grid,franke.fn,franke.data)
importFrom("sp","coordinates")
importFrom("sp","coordinates<-")
......
R/interp.R
View file @ b9de5f5f
......@@ -25,6 +25,11 @@ interp <-
} else
stop("either x,y,z are numerical or x is SpatialPointsDataFrame and z a name of a data column in x")
}
## FIXME: drop old akima code
if(linear)
ret <- interp.old(x,y,z,xo,yo,ncp=0,extrap=FALSE,
duplicate=duplicate,dupfun=dupfun)
else {
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
......@@ -68,8 +73,8 @@ interp <-
n <- length(x)
}
zo <- matrix(0, nx, ny)
storage.mode(zo) <- "double"
## zo <- matrix(0, nx, ny)
## storage.mode(zo) <- "double"
miss <- !extrap # if not extrapolating, set missing values
ans <- .Fortran("sdsf3p",
......@@ -82,7 +87,7 @@ interp <-
x = as.double(xo),
ny=as.integer(ny),
y = as.double(yo),
z = zo,
z = as.double(matrix(0,nx,ny)),
ier = integer(1),
wk = double(36 * n),
iwk = integer(25 * n),
......@@ -143,7 +148,7 @@ interp <-
x = as.double(xo),
as.integer(ny),
y = as.double(yo),
z = zo,
z = as.double(matrix(0,nx,ny)),
ier = integer(1),
double(36 * n),
integer(25 * n),
......@@ -177,7 +182,8 @@ interp <-
ret@proj4string <- sp.proj4string
gridded(ret) <- TRUE
} else {
ret <- list(x=ans$x,y=ans$y,z=ans$z)
ret <- list(x=ans$x,y=ans$y,z=matrix(ans$z,nx,ny))
}
} ## END FIXME
ret
}
R/interp.new.R 0 → 100644
View file @ b9de5f5f
interp.new <-
function(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)
{
warning("interp.new() is deprecated, use interp()")
interp(x, y, z,
xo,
yo, linear,
extrap, duplicate, dupfun)
}
R/interp.old.R 0 → 100644
View file @ b9de5f5f
"interp.old"<-function(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)
{
## warning("interp.old() is deprecated, future versions will only provide interp()")
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
if(ncp > 25) {
ncp <- 25
cat("ncp too large, using ncp=25\n")
}
drx <- diff(range(x))
dry <- diff(range(y))
if(drx == 0 || dry == 0)
stop("all data collinear") # other cases caught in Fortran code
if(drx/dry > 10000 || drx/dry < 0.0001)
stop("scales of x and y are too dissimilar")
n <- length(x)
nx <- length(xo)
ny <- length(yo)
if(length(y) != n || length(z) != n)
stop("Lengths of x, y, and z do not match")
xy <- paste(x, y, sep = ",")# trick for 'duplicated' (x,y)-pairs
if(duplicate == "error") {
if(any(duplicated(xy)))
stop("duplicate data points: need to set 'duplicate = ..' ")
}
else { ## duplicate != "error"
i <- match(xy, xy)
if(duplicate == "user")
dupfun <- match.fun(dupfun)#> error if it fails
ord <- !duplicated(xy)
if(duplicate != "strip") {
centre <- function(x)
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x))
z <- unlist(lapply(split(z,i), centre))
} else {
z <- z[ord]
}
x <- x[ord]
y <- y[ord]
n <- length(x)
}
zo <- matrix(0, nx, ny)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
misso <- matrix(miss, nx, ny)
if(extrap & ncp == 0)
warning("Cannot extrapolate with linear option")
ans <- .Fortran("idsfft",
as.integer(1),
as.integer(ncp),
as.integer(n),
as.double(x),
as.double(y),
as.double(z),
as.integer(nx),
as.integer(ny),
x = as.double(xo),
y = as.double(yo),
z = zo,
integer((31 + ncp) * n + nx * ny),
double(5 * n),
misso = as.logical(misso),
PACKAGE = "akima")[c("x", "y", "z", "misso")]
ans$z[ans$misso] <- NA
ans[c("x", "y", "z")]
}
R/interpp.R
View file @ b9de5f5f
......@@ -29,6 +29,12 @@
"or both x and xo have to be SpatialPointsDataFrames",
"and z a name of a data column in x"))
}
## FIXME: drop old akima code
if(linear)
ret <- interpp.old(x,y,z,xo,yo,ncp=0,extrap=FALSE,
duplicate=duplicate,dupfun=dupfun)
else {
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
if(is.null(xo))
......@@ -164,5 +170,6 @@
} else {
ret <- list(x=ans$x,y=ans$y,z=ans$z)
}
} ## END FIXME
ret
}
R/interpp.new.R 0 → 100644
View file @ b9de5f5f
"interpp.new"<-function(x, y, z, xo, yo, extrap = FALSE,
duplicate = "error", dupfun = NULL)
{
warning("interp.new() is deprecated, use interp()")
interpp(x, y, z,
xo,
yo, linear=FALSE,
extrap=extrap, duplicate=duplicate, dupfun=dupfun)
}
R/interpp.old.R 0 → 100644
View file @ b9de5f5f
"interpp.old"<-function(x, y, z, xo, yo, ncp = 0, extrap = FALSE,
duplicate = "error", dupfun = NULL)
{
warning("interpp.old() is deprecated, future versions will only provide interpp()")
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
if(is.null(xo))
stop("xo missing")
if(is.null(yo))
stop("yo missing")
if(ncp>25){
ncp <- 25
cat("ncp too large, using ncp=25\n")
}
drx <- diff(range(x))
dry <- diff(range(y))
if(drx == 0 || dry == 0)
stop("all data collinear") # other cases caught in Fortran code
if(drx/dry > 10000 || drx/dry < 0.0001)
stop("scales of x and y are too dissimilar")
n <- length(x)
np <- length(xo)
if(length(yo)!=np)
stop("length of xo and yo differ")
if(length(y) != n || length(z) != n)
stop("Lengths of x, y, and z do not match")
xy <- paste(x, y, sep =",")
i <- match(xy, xy)
if(duplicate=="user" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x)
)
}
if(duplicate!="strip"){
z <- unlist(lapply(split(z,i), centre))
ord <- !duplicated(xy)
x <- x[ord]
y <- y[ord]
n <- length(x)
}
else{
ord <- (hist(i,plot=FALSE,freq=TRUE,breaks=seq(0.5,max(i)+0.5,1))$counts==1)
x <- x[ord]
y <- y[ord]
z <- z[ord]
n <- length(x)
}
}
else
if(any(duplicated(xy)))
stop("duplicate data points")
zo <- rep(0, np)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
misso <- seq(miss, np)
if(extrap & ncp == 0)
warning("Cannot extrapolate with linear option")
ans <- .Fortran("idbvip",
as.integer(1),
as.integer(ncp),
as.integer(n),
as.double(x),
as.double(y),
as.double(z),
as.integer(np),
x = as.double(xo),
y = as.double(yo),
z = zo,
integer((31 + ncp) * n + np),
double(8 * n),
misso = as.logical(misso),
PACKAGE = "akima")
temp <- ans[c("x", "y", "z", "misso")]
temp$z[temp$misso]<-NA
temp[c("x", "y", "z")]
}
man/interp.old.Rd 0 → 100644
View file @ b9de5f5f
\name{interp.old}
\title{Gridded Bivariate Interpolation for Irregular Data}
\alias{interp.new}
\alias{interp.old}
\description{
These functions implement bivariate interpolation onto a grid
for irregularly spaced input data. These functions are only for
backward compatibility, use \code{\link{interp}} instead.
}
\usage{
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{
\item{x}{
vector of x-coordinates of data points or a
\code{SpatialPointsDataFrame} object.
Missing values are not accepted.
}
\item{y}{
vector of y-coordinates of data points.
Missing values are not accepted.
If left as NULL indicates that \code{x} should be a
\code{SpatialPointsDataFrame} and \code{z} names the variable of
interest in this dataframe.
}
\item{z}{
vector of z-coordinates of data points or a character variable
naming the variable of interest in the
\code{SpatialPointsDataFrame} \code{x}.
Missing values are not accepted.
\code{x}, \code{y}, and \code{z} must be the same length
(execpt if \code{x} is a \code{SpatialPointsDataFrame}) and may
contain no fewer than four points. The points of \code{x} and
\code{y} cannot be collinear, i.e, they cannot fall on the same line
(two vectors \code{x} and \code{y} such that \code{y = ax + b} for
some \code{a}, \code{b} will not be accepted). \code{interp} is
meant for cases in which you have \code{x}, \code{y} values
scattered over a plane and a \code{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 \code{outer()}.
}
\item{xo}{
vector of x-coordinates of output grid. The default is 40 points
evenly spaced over the range of \code{x}. If extrapolation is not being
used (\code{extrap=FALSE}, the default), \code{xo} should have a
range that is close to or inside of the range of \code{x} for the
results to be meaningful.
}
\item{yo}{vector of y-coordinates of output grid; analogous to
\code{xo}, see above.}
\item{linear}{logical -- indicating wether linear or spline
interpolation should be used. supersedes old \code{ncp} parameter}
\item{ncp}{
deprecated, use parameter \code{linear}. Now only used by
\code{interp.old()}.
old meaning was:
number of additional points to be used in computing partial
derivatives at each data point.
\code{ncp} must be either \code{0} (partial derivatives are not used), or at
least 2 but smaller than the number of data points (and smaller than
25).
}
\item{extrap}{
logical flag: should extrapolation be used outside of the
convex hull determined by the data points?}
\item{duplicate}{character string indicating how to handle duplicate
data points. Possible values are
\describe{
\item{\code{"error"}}{produces an error message,}
\item{\code{"strip"}}{remove duplicate z values,}
\item{ \code{"mean"},\code{"median"},\code{"user"}}{calculate
mean , median or user defined function (\code{dupfun}) of duplicate
z values.}
}}
\item{dupfun}{a function, applied to duplicate points if
\code{duplicate= "user"}.}
}
\value{
list with 3 components:
\item{x,y}{
vectors of x- and y- coordinates of output grid, the same as the input
argument \code{xo}, or \code{yo}, if present. Otherwise, their
default, a vector 40 points evenly spaced over the range of the
input \code{x}.}
\item{z}{
matrix of fitted z-values. The value \code{z[i,j]} is computed
at the x,y point \code{xo[i], yo[j]}. \code{z} has
dimensions \code{length(xo)} times \code{length(yo)}.}
If input is a \code{SpatialPointsDataFrame} a
\code{SpatialPixelssDataFrame} is returned.
}
\note{
\code{interp.new} is deprecated and \code{interp.old} will soon be
deprecated.
}
\details{
see \code{\link{interp}}
}
\references{
Akima, H. (1978). A Method of Bivariate Interpolation and
Smooth Surface Fitting for Irregularly Distributed Data Points.
ACM Transactions on Mathematical Software \bold{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 \bold{22}, 362--371.
}
\seealso{
\code{\link{contour}}, \code{\link{image}},
\code{\link{approx}}, \code{\link{spline}},
\code{\link{aspline}},
\code{\link{outer}}, \code{\link{expand.grid}}.
}
\keyword{dplot}
man/interpp.Rd
View file @ b9de5f5f
......@@ -105,6 +105,10 @@ interpp(x, y=NULL, z, xo, yo=NULL, linear=TRUE, extrap=FALSE,
See \code{\link{interp}} for more details.
}
\description{
These functions implement bivariate interpolation onto a set of points
for irregularly spaced input data. These functions are only for
backward compatibility, use \code{\link{interpp}} instead.
If \code{linear} is \\code{TRUE}, linear
interpolation is used in the triangles bounded by data points, otherwise
cubic interpolation is done.
......
man/interpp.old.Rd 0 → 100644
View file @ b9de5f5f
\name{interpp.old}
\title{
Pointwise Bivariate Interpolation for Irregular Data
}
\alias{interpp.old}
\alias{interpp.new}
\usage{
interpp.old(x, y, z, xo, yo, ncp = 0, extrap = FALSE,
duplicate = "error", dupfun = NULL)
interpp.new(x, y, z, xo, yo, extrap = FALSE,
duplicate = "error", dupfun = NULL)
}
\arguments{
\item{x}{
vector of x-coordinates of data points or a
\code{SpatialPointsDataFrame} object.
Missing values are not accepted.
}
\item{y}{
vector of y-coordinates of data points.
Missing values are not accepted.
If left as NULL indicates that \code{x} should be a
\code{SpatialPointsDataFrame} and \code{z} names the variable of
interest in this dataframe.
}
\item{z}{
vector of z-coordinates of data points or a character variable
naming the variable of interest in the
\code{SpatialPointsDataFrame} \code{x}.
Missing values are not accepted.
\code{x}, \code{y}, and \code{z} must be the same length
(execpt if \code{x} is a \code{SpatialPointsDataFrame}) and may contain no fewer
than four points. The points of \code{x} and \code{y}
cannot be collinear, i.e, they cannot fall on the same line (two vectors
\code{x} and \code{y} such that \code{y = ax + b} for some \code{a}, \code{b} will not be
accepted).
}
\item{xo}{
vector of x-coordinates of points at which to evaluate the interpolating
function. If \code{x} is a \code{SpatialPointsDataFrame} this has
also to be a \code{SpatialPointsDataFrame}.
}
\item{yo}{
vector of y-coordinates of points at which to evaluate the interpolating
function.
If operating on \code{SpatialPointsDataFrame}s this is left as \code{NULL}
}
\item{ncp}{
deprecated, use parameter \code{linear}. Now only used by
\code{interpp.old()}.
meaning was:
number of additional points to be used in computing partial
derivatives at each data point.
\code{ncp} must be either \code{0} (partial derivatives are not used, =
linear interpolation), or at
least 2 but smaller than the number of data points (and smaller than 25).
}
\item{extrap}{
logical flag: should extrapolation be used outside of the
convex hull determined by the data points?}
\item{duplicate}{
indicates how to handle duplicate data points. Possible values are
\code{"error"} - produces an error message, \code{"strip"} - remove
duplicate z values, \code{"mean"},\code{"median"},\code{"user"} -
calculate mean , median or user defined function of duplicate z
values.
}
\item{dupfun}{this function is applied to duplicate points if \code{duplicate="user"}
}
}
\value{
list with 3 components:
\item{x}{
vector of x-coordinates of output points, the same as the input
argument \code{xo}.
}
\item{y}{
vector of y-coordinates of output points, the same as the input
argument \code{yo}.
}
\item{z}{
fitted z-values. The value \code{z[i]} is computed
at the x,y point \code{x[i], y[i]}.
}
If input is \code{SpatialPointsDataFrame} than an according
\code{SpatialPointsDataFrame} is returned.
}
\section{NOTE}{
Use \code{interp} if interpolation on a regular grid is wanted.
The two versions \code{interpp.old} and \code{interpp.new} are now
deprecated, use \code{\link{interpp}} instead, see details there.
Earlier versions (pre 0.5-1) of \code{interpp} used the parameter
\code{ncp} to choose between linear and cubic interpolation, this is now done
by setting the logical parameter \code{linear}. Use of \code{ncp} is still
possible, but is deprecated.
}
\description{
If \code{ncp} is zero, linear
interpolation is used in the triangles bounded by data points.
Cubic interpolation is done if partial derivatives are used.
If \code{extrap} is \code{FALSE}, z-values for points outside the convex hull are
returned as \code{NA}.
No extrapolation can be performed if \code{ncp} is zero.
The \code{interpp} function handles duplicate \code{(x,y)} points in
different ways. As default it will stop with an error message. But it
can give duplicate points an unique \code{z} value according to the
parameter \code{duplicate} (\code{mean},\code{median} or any other
user defined function).
The triangulation scheme used by \code{interp} works well if \code{x} and \code{y} have
similar scales but will appear stretched if they have very different
scales. The spreads of \code{x} and \code{y} must be within four orders of magnitude
of each other for \code{interpp} 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,
\bold{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,
\bold{22}, 362-371.
}
\seealso{
\code{\link[graphics]{contour}}, \code{\link[graphics]{image}},
\code{\link[stats]{approxfun}}, \code{\link[stats]{splinefun}},
\code{\link[base]{outer}}, \code{\link[base]{expand.grid}},
\code{\link{interp}}, \code{\link{aspline}}.
}
\keyword{dplot}
src/akima.h
View file @ b9de5f5f
#include <R.h>
/* ACM 526, soon to be removed: */
int F77_NAME(idbvip) (int *md, int *ncp, int *ndp,
double *xd, double *yd, double *zd,
int *nip, double *xi, double *yi, double *zi,
int *iwk, double *wk, int *missi);
int F77_NAME(idsfft) (int *md, int *ncp, int *ndp,
double *xd, double *yd,double *zd,
int *nxi, int *nyi,
double *xi, double *yi, double *zi,
int *iwk, double *wk, int *missi);
/* ACM 679: */
int F77_NAME(uvip3p) (int *np, int *nd, double *xd, double *yd,
......@@ -25,6 +36,6 @@ int F77_NAME(sdbi3p) (int *md, int *ndp, double *xd, double *yd, double *zd,
int *extrpi, int *near, int *next, double *dist);
/* bilinear, A. Gebhardt: */
int F77_NAME(biliip) (double *x0, double *y0, double *z0,
int F77_NAME(biliip) (double *x0, double *y0, double *z0,
double *x, double *y, double *z,
int *nx, int *ny);
src/idbvip.f 0 → 100644
View file @ b9de5f5f
SUBROUTINE IDBVIP(MD,NCP,NDP,XD,YD,ZD,NIP,XI,YI,ZI, ID001340
1 IWK,WK,MISSI)
C THIS SUBROUTINE PERFORMS BIVARIATE INTERPOLATION WHEN THE PRO-
C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY
C DISTRIBUTED IN THE PLANE.
C THE INPUT PARAMETERS ARE
C MD = MODE OF COMPUTATION (MUST BE 1, 2, OR 3),
C = 1 FOR NEW NCP AND/OR NEW XD-YD,
C = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI,
C = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI,
C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C MATING PARTIAL DERIVATIVES AT EACH DATA POINT
C (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP),
C NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER),
C XD = ARRAY OF DIMENSION NDP CONTAINING THE X
C COORDINATES OF THE DATA POINTS,
C YD = ARRAY OF DIMENSION NDP CONTAINING THE Y