Skip to content

GitLab

Commit b7940cf6 authored by agebhard's avatar agebhard
Browse files

bugfixes/enhancements from Martin Maechler (thanks :-)

parent c1322d6b
Show whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ b7940cf6
Package: akima
Version: 0.4-4
Version: 0.4-5
Date: 2005-07-21
Title: Interpolation of irregularly spaced data
Author: Fortran code by H. Akima
R port by Albrecht Gebhardt <albrecht.gebhardt@uni-klu.ac.at>
aspline function by Thomas Petzoldt <petzoldt@rcs.urz.tu-dresden.de>
enhancements and corrections by Martin Maechler <maechler@stat.math.ethz.ch>
Maintainer: Albrecht Gebhardt <albrecht.gebhardt@uni-klu.ac.at>
Description: Linear or cubic spline interpolation for irregular gridded data
License: Fortran code: ACM, free for non-commercial use, R functions GPL
R/interp.R
View file @ b7940cf6
"interp"<-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)
interp <-
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)
{
if (ncp==0)
# use the old version for linear interpolation
interp.old(x, y, z, xo, yo, ncp, extrap, duplicate, dupfun)
else
interp.new(x, y, z, xo, yo, ncp, extrap, duplicate, dupfun)
if (ncp == 0)
## use the old version for linear interpolation
interp.old(x, y, z, xo = xo, yo = yo, ncp = ncp,
extrap = extrap, duplicate = duplicate, dupfun = dupfun)
else ## use the new one
interp.new(x, y, z, xo = xo, yo = yo, linear = FALSE,
ncp = NULL,# not using 'ncp' argument
extrap = extrap, duplicate = duplicate, dupfun = dupfun)
}
R/interp.new.R
View file @ b7940cf6
"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)
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)
{
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
if(!is.null(ncp)){
if(ncp!=0){
if(!is.null(ncp)) {
if(ncp != 0) {
cat("ncp not supported, it is automatically choosen by Fortran code\n")
}
else {
......@@ -14,7 +15,7 @@
stop("use interp.old().")
}
}
if(linear){
if(linear) {
cat("linear interpolation not yet implemented with interp.new().\n")
stop("use interp.old().")
}
......@@ -30,48 +31,42 @@
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 =",")
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" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
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))
}
if(duplicate!="strip"){
z <- unlist(lapply(split(z,i), centre))
ord <- !duplicated(xy)
x <- x[ord]
y <- y[ord]
n <- length(x)
} else {
z <- z[ord]
}
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 <- matrix(0, nx, ny)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
extrap <- matrix(TRUE, nx, ny)
if(!is.null(ncp)){
if(extrap & ncp == 0)
warning("Cannot extrapolate with linear option")
}
else {
if(extrap & linear)
miss <- !extrap # if not extrapolating, set missing values
misso <- matrix(TRUE, nx, ny)# hmm, or rather 'miss' ??
if(extrap && if(is.null(ncp)) linear else (ncp == 0))
warning("Cannot extrapolate with linear option")
}
ans <- .Fortran("sdsf3p",
as.integer(1),
as.integer(n),
......@@ -86,13 +81,13 @@
ier = integer(1),
double(36 * n),
integer(25 * n),
extrap = as.logical(extrap),
extrap = as.logical(misso),
near = integer(n),
nxt = integer(n),
dist = double(n),
PACKAGE = "akima")
temp <- ans[c("x", "y", "z", "extrap")]
PACKAGE = "akima")[c("x", "y", "z", "extrap")]
if(miss)
temp$z[temp$extrap]<-NA
temp[c("x", "y", "z")]
ans$z[ans$extrap] <- NA
ans[c("x", "y", "z")]
}
R/interp.old.R
View file @ b7940cf6
......@@ -4,10 +4,11 @@
{
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
if(ncp>25){
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)
......@@ -19,36 +20,34 @@
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 =",")
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" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
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))
}
if(duplicate!="strip"){
z <- unlist(lapply(split(z,i), centre))
ord <- !duplicated(xy)
x <- x[ord]
y <- y[ord]
n <- length(x)
} else {
z <- z[ord]
}
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 <- matrix(0, nx, ny)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
......@@ -70,8 +69,7 @@
integer((31 + ncp) * n + nx * ny),
double(5 * n),
misso = as.logical(misso),
PACKAGE = "akima")
temp <- ans[c("x", "y", "z", "misso")]
temp$z[temp$misso]<-NA
temp[c("x", "y", "z")]
PACKAGE = "akima")[c("x", "y", "z", "misso")]
ans$z[ans$misso] <- NA
ans[c("x", "y", "z")]
}
man/interp.Rd
View file @ b7940cf6
\name{interp}
\title{
Gridded Bivariate Interpolation for Irregular Data
}
\title{Gridded Bivariate Interpolation for Irregular Data}
\alias{interp}
\alias{interp.new}
\alias{interp.old}
\description{
These functions implement bivariate interpolation onto a grid
for irregularly spaced input data. Bilinear or bicubic spline
interpolation is applied using different versions of algorithms from
Akima.
}
\usage{
interp(x, y, z, xo=seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), length = 40), ncp=0, extrap=FALSE, duplicate = "error", dupfun = NULL)
interp.old(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), length = 40), ncp=0, extrap=FALSE, duplicate = "error", dupfun = NULL)
interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), length = 40), linear=FALSE, ncp=NULL, extrap=FALSE, duplicate = "error", dupfun = NULL)
interp(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.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}{
......@@ -23,30 +33,27 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y),
vector of z-coordinates of data points.
Missing values are not accepted.
\code{x}, \code{y}, and \code{z} must be the same length 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()}.
\code{x}, \code{y}, and \code{z} must be the same length 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. The default is 40 points
evenly spaced over the range of \code{y}. If extrapolation is not being
used (\code{extrap=FALSE}, the default), \code{yo} should have a range that is
close to or inside of the range of \code{y} for the results to be meaningful.
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{linear}{logical, switch to linear interpolation in \code{interp.new}}
\item{yo}{vector of y-coordinates of output grid; analogous to
\code{xo}, see above.}
\item{linear}{logical, switch to linear interpolation in \code{interp.new}.}
\item{ncp}{
number of additional points to be used in computing partial
derivatives at each data point.
......@@ -56,38 +63,31 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y),
}
\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"}
}
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}{
vector of x-coordinates of output grid, the same as the input
argument \code{xo}, if present. Otherwise, a vector 40 points evenly spaced
over the range of the input \code{x}.
}
\item{y}{
vector of y-coordinates of output grid, the same as the input
argument \code{yo}, if present. Otherwise, a vector 40 points evenly spaced
over the range of the input \code{x}.
}
\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{x[i], y[j]}. \code{z} has
dimensions \code{length(x)} times \code{length(y)} (\code{length(xo)} times \code{length(yo)}).
}}
at the x,y point \code{xo[i], yo[j]}. \code{z} has
dimensions \code{length(xo)} times \code{length(yo)}.}
}
\note{
\code{interp} is a wrapper for the two versions \code{interp.old} (it
uses (almost) the same Fortran code from Akima 1978 as the S-Plus version) and
......@@ -103,7 +103,7 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y),
functions \code{contour} and \code{image}. Check the requirements of
these functions when choosing values for \code{xo} and \code{yo}.
}
\description{
\details{
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.
......@@ -111,7 +111,6 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y),
convex hull are returned as \code{NA}.
No extrapolation can be performed if \code{ncp} is zero.
The \code{interp} 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
......@@ -137,50 +136,75 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y),
}
\seealso{
\code{\link[base]{contour}}, \code{\link[base]{image}}, \code{\link[base]{approx}}, \code{\link[base]{spline}}, \code{\link[base]{outer}}, \code{\link[base]{expand.grid}}.
\code{\link{contour}}, \code{\link{image}},
\code{\link{approx}}, \code{\link{spline}},
\code{\link{outer}}, \code{\link{expand.grid}}.
}
\examples{
data(akima)
# linear interpolation
plot(y ~ x, data = akima, main = "akima example data")
with(akima, text(x, y, formatC(z,dig=2), adj = -0.1))
## linear interpolation
akima.li <- interp(akima$x, akima$y, akima$z)
image(akima.li$x,akima.li$y,akima.li$z)
contour(akima.li$x,akima.li$y,akima.li$z,add=TRUE)
points(akima$x,akima$y)
# increase smoothness
akima.smooth <- interp(akima$x, akima$y, akima$z,
xo=seq(0,25, length=100), yo=seq(0,20, length=100))
image(akima.smooth$x,akima.smooth$y,akima.smooth$z)
contour(akima.smooth$x,akima.smooth$y,akima.smooth$z,add=TRUE)
points(akima$x,akima$y)
# use triangulation library to
# show underlying triangulation:
image (akima.li, add=TRUE)
contour(akima.li, add=TRUE)
points (akima, pch = 3)
## increase smoothness (using finer grid):
akima.smooth <-
with(akima, interp(x, y, z, xo=seq(0,25, length=100),
yo=seq(0,20, length=100)))
image (akima.smooth)
contour(akima.smooth, add=TRUE)
points(akima, pch = 3, cex = 2, col = "blue")
# use triangulation package to show underlying triangulation:
if(library(tripack, logical.return=TRUE))
plot(tri.mesh(akima),add=TRUE,lty="dashed")
plot(tri.mesh(akima), add=TRUE, lty="dashed")
# use only 15 points (interpolation only within convex hull!)
akima.part <- interp(akima$x[1:15],akima$y[1:15],akima$z[1:15])
image(akima.part$x,akima.part$y,akima.part$z)
contour(akima.part$x,akima.part$y,akima.part$z,add=TRUE)
akima.part <- with(akima, interp(x[1:15], y[1:15], z[1:15]))
image (akima.part)
contour(akima.part, add=TRUE)
points(akima$x[1:15],akima$y[1:15])
# spline interpolation, use 5 points to calculate derivatives
# interp gives `linear interpolation not yet implemented with interp.new()'
akima.spl <- interp.old(akima$x, akima$y, akima$z,
xo=seq(0,25, length=100), yo=seq(0,20, length=100),ncp=5)
image(akima.spl$x,akima.spl$y,akima.spl$z)
contour(akima.spl$x,akima.spl$y,akima.spl$z,add=TRUE)
points(akima$x,akima$y)
## spline interpolation
## --------------------
## "Old": use 5 points to calculate derivatives -> many NAs
akima.sO <- interp.old(akima$x, akima$y, akima$z,
xo=seq(0,25, length=100), yo=seq(0,20, length=100), ncp=5)
table(is.na(akima.sO$z)) ## 3990 NA's; = 40 \%
akima.sO <- with(akima,
interp.old(x,y,z, xo=seq(0,25, length=100), yo=seq(0,20, len=100), ncp = 4))
sum(is.na(akima.sO$z)) ## still 3429
image (akima.sO) # almost useless
contour(akima.sO, add = TRUE)
## "New:"
akima.spl <- with(akima, interp.new(x,y,z, xo=seq(0,25, length=100),
yo=seq(0,20, length=100)))
contour(akima.spl) ; points(akima)
full.pal <- function(n) hcl(h = seq(340, 20, length = n))
cool.pal <- function(n) hcl(h = seq(120, 0, length = n) + 150)
warm.pal <- function(n) hcl(h = seq(120, 0, length = n) - 30)
filled.contour(akima.spl, color.palette = full.pal,
plot.axes = { axis(1); axis(2);
points(akima, pch = 3, col= hcl(c=100, l = 20))})
# no extrapolation!
## example with duplicate points :
# example with duplicate points
data(airquality)
air <- airquality[(!is.na(airquality$Temp) &
!is.na(airquality$Ozone) &
!is.na(airquality$Solar.R)),]
# gives an error:
\dontrun{air.ip <- interp(air$Temp,air$Solar.R,air$Ozone)}
air <- subset(airquality,
!is.na(Temp) & !is.na(Ozone) & !is.na(Solar.R))
# gives an error {duplicate ..}:
try( air.ip <- interp.new(air$Temp,air$Solar.R,air$Ozone) )
# use mean of duplicate points:
air.ip <- interp(air$Temp,air$Solar.R,air$Ozone,duplicate="mean")
air.ip <- with(air, interp.new(Temp, Solar.R, log(Ozone), duplicate = "mean"))
image(air.ip, main = "Airquality: Ozone vs. Temp and Solar.R")
with(air, points(Temp, Solar.R))
}
\keyword{dplot}
% Converted by Sd2Rd version 0.2-a3.
orig/akima.bib
View file @ b7940cf6
<HTML>
<TITLE>TOMS : Bibliographic record for "Akima:1996:ASS" </TITLE>
<BODY BGCOLOR="#fffaf0">
<H1><A HREF="/toms/"><img align=bottom src="/toms/toms.gif"></A>
Bibliographic record for "Akima:1996:ASS"</H1>
<P></P><HR><P></P>
<PRE>
@Article{Akima:1996:ASS,
author = "Hiroshi Akima",
title = "Algorithm 761: scattered-data surface fitting that has
......@@ -30,8 +23,3 @@ Bibliographic record for "Akima:1996:ASS"</H1>
G.4}: Mathematics of Computing, MATHEMATICAL
SOFTWARE.",
}
</PRE>
<P></P><HR><P></P>
</BODY>
</HTML>
src/idptip.f
View file @ b7940cf6
......@@ -3,6 +3,7 @@
1 ZII,MISSII)
C THIS SUBROUTINE PERFORMS PUNCTUAL INTERPOLATION OR EXTRAPOLA-
C TION, I.E., DETERMINES THE Z VALUE AT A POINT.
C THE INPUT PARAMETERS ARE
C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X,
C Y, AND Z COORDINATES OF THE DATA POINTS, WHERE
......@@ -22,9 +23,11 @@ C THE POINT FOR WHICH INTERPOLATION IS TO BE
C PERFORMED,
C XII,YII = X AND Y COORDINATES OF THE POINT FOR WHICH
C INTERPOLATION IS TO BE PERFORMED.
C THE OUTPUT PARAMETER IS
C THE OUTPUT PARAMETERs are
C ZII = INTERPOLATED Z VALUE.
C MISSII = LOCIGAL INDICATING MISSING VALUE
C DECLARATION STATEMENTS
IMPLICIT DOUBLE PRECISION (A-D,P-Z)
LOGICAL MISSII
......@@ -35,6 +38,7 @@ C DECLARATION STATEMENTS
1 ZU(3),ZV(3),ZUU(3),ZUV(3),ZVV(3)
DOUBLE PRECISION LU,LV
EQUIVALENCE (P5,P50)
C PRELIMINARY PROCESSING
10 IT0=ITI
NTL=NT+NL
......@@ -170,6 +174,7 @@ C EVALUATES THE POLYNOMIAL.
ZII=P0+U*(P1+U*(P2+U*(P3+U*(P4+U*P5))))
MISSII=.FALSE.
RETURN
C CALCULATION OF ZII BY EXTRAPOLATION IN THE RECTANGLE.
C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED.
40 IF(IT0.EQ.ITPV) GO TO 50
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment