Commit 490a10df authored by agebhard's avatar agebhard
Browse files

* Fixes for R 3.4.0 (gfortran specifc)

* old Akima code removed
* interp.old removed
* interp and interp.new merged
* similar for interpp
* a bilinear interpolation fro rectangular grids added (also gridwise in bilinear.grid)
* new ACM license url
parent 6021ca47
Package: akima
Version: 0.5-12
Date: 2015-08-26
Version: 0.6-1
Date: 2016-12-09
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)"),
......
This package is distributed under the ACM license at
http://www.acm.org/publications/policies/softwarecrnotice.
http://www.acm.org/publications/policies/software-copyright-notice
In detail:
1. Fortran code (src/*.f):
Copyrighted and Licensed by ACM,
see http://www.acm.org/publications/policies/softwarecrnotice
see http://www.acm.org/publications/policies/software-copyright-notice
2. R interface (src/init.c src/akima.h R/*.R man/*.Rd data/*):
......
useDynLib(akima)
export(aspline,interp,interpp,interp.old,interp.new,interp2xyz,bicubic,bicubic.grid)
export(aspline,interp,interpp,interp2xyz,bicubic,bicubic.grid,bilinear,bilinear.grid)
importFrom("sp","coordinates")
importFrom("sp","coordinates<-")
......@@ -9,3 +9,4 @@ importFrom("sp","gridded<-")
importFrom("grDevices", "xy.coords")
importFrom("graphics", "hist")
importFrom("stats", "median")
importFrom("stats", "rnorm")
......@@ -31,19 +31,27 @@ bicubic <- function(x,y,z,x0,y0){
list(x=x0,y=y0,z=ret$zi)
}
bicubic.grid <- function(x,y,z,xlim,ylim,dx,dy){
nx <- length(x)
ny <- length(y)
if(dim(z)[1]!=nx)
bicubic.grid <- function(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)),
nx=40,ny=40,dx=NULL,dy=NULL){
Nx <- length(x)
Ny <- length(y)
if(dim(z)[1]!=Nx)
stop("dim(z)[1] and length of x differs!")
if(dim(z)[2]!=ny)
if(dim(z)[2]!=Ny)
stop("dim(z)[2] and length of y differs!")
if(!is.null(dx)){
xi <- seq(xlim[1],xlim[2],by=dx)
nx <- length(xi)
} else {
xi <- seq(ylim[1],ylim[2],length=nx)
}
xi <- seq(xlim[1],xlim[2],by=dx)
nx <- length(xi)
yi <- seq(ylim[1],ylim[2],by=dy)
ny <- length(yi)
if(!is.null(dx)){
yi <- seq(ylim[1],ylim[2],by=dy)
ny <- length(yi)
} else {
yi <- seq(ylim[1],ylim[2],length=ny)
}
xmat <- matrix(rep(xi,ny),nrow=ny,ncol=nx,byrow=TRUE)
ymat <- matrix(rep(yi,nx),nrow=ny,ncol=nx,byrow=FALSE)
......
bilinear <- function(x,y,z,x0,y0){
nx <- length(x)
ny <- length(y)
if(dim(z)[1]!=nx)
stop("dim(z)[1] and length of x differs!")
if(dim(z)[2]!=ny)
stop("dim(z)[2] and length of y differs!")
n0 <- length(x0)
if(length(y0)!=n0)
stop("length of y0 and x0 differs!")
ret <- .Fortran("biliip",
as.double(x0),
as.double(y0),
z0=double(n0),
as.integer(n0),
as.double(x),
as.double(y),
as.double(z),
as.integer(nx),
as.integer(ny),
PACKAGE="akima")
list(x=x0,y=y0,z=ret$z0)
}
bilinear.grid <- function(x,y,z,xlim=c(min(x),max(x)),ylim=c(min(y),max(y)),
nx=40,ny=40,dx=NULL,dy=NULL){
Nx <- length(x)
Ny <- length(y)
if(dim(z)[1]!=Nx)
stop("dim(z)[1] and length of x differs!")
if(dim(z)[2]!=Ny)
stop("dim(z)[2] and length of y differs!")
if(!is.null(dx)){
xi <- seq(xlim[1],xlim[2],by=dx)
nx <- length(xi)
} else {
xi <- seq(ylim[1],ylim[2],length=nx)
}
if(!is.null(dx)){
yi <- seq(ylim[1],ylim[2],by=dy)
ny <- length(yi)
} else {
yi <- seq(ylim[1],ylim[2],length=ny)
}
xmat <- matrix(rep(xi,ny),nrow=ny,ncol=nx,byrow=TRUE)
ymat <- matrix(rep(yi,nx),nrow=ny,ncol=nx,byrow=FALSE)
xy <- cbind(c(xmat),c(ymat))
n0 <- nx*ny
ret <- bilinear(x,y,z,xy[,1],xy[,2])
# return cell boundaries
list(x=xi,y=yi,z=t(matrix(ret$z,nrow=ny,ncol=nx,byrow=F)))
}
interp <-
function(x, y=NULL, z,
xo = seq(min(x), max(x), length = nx),
yo = seq(min(y), max(y), length = ny), linear=TRUE,
extrap = FALSE, duplicate = "error", dupfun = NULL, ncp=NULL,
nx=40, ny=40)
function(x, y=NULL, z,
xo = seq(min(x), max(x), length = nx),
yo = seq(min(y), max(y), length = ny), linear = TRUE,
extrap = FALSE, duplicate = "error", dupfun = NULL,
nx=40, ny=40,
jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE)
{
## for backward compatibility
if(!is.null(ncp)){
warning('use of \'ncp\' parameter is deprecated!')
if(ncp==0)
linear <- TRUE
else if(ncp>0)
linear <- FALSE
else
stop('ncp < 0 ?')
}
## handle sp data, save coordinate and value names
is.sp <- FALSE
sp.coord <- NULL
......@@ -36,29 +26,137 @@ interp <-
stop("either x,y,z are numerical or x is SpatialPointsDataFrame and z a name of a data column in x")
}
if(linear)
## use the old version for linear interpolation
ret <- interp.old(x, y, z, xo = xo, yo = yo, ncp = 0,
extrap = extrap, duplicate = duplicate,
dupfun = dupfun)
else ## use the new one
ret <- interp.new(x, y, z, xo = xo, yo = yo, linear = FALSE,
ncp = NULL,# not using 'ncp' argument
extrap = extrap, duplicate = duplicate,
dupfun = dupfun)
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
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, set missing values
ans <- .Fortran("sdsf3p",
md = as.integer(1),
ndp = as.integer(n),
xd = as.double(x),
yd = as.double(y),
zd = as.double(z),
nx = as.integer(nx),
x = as.double(xo),
ny=as.integer(ny),
y = as.double(yo),
z = zo,
ier = integer(1),
wk = double(36 * n),
iwk = integer(25 * n),
extrap = as.logical(matrix(extrap,nx,ny)),
near = integer(n),
nxt = integer(n),
dist = double(n),
linear = as.logical(linear),
PACKAGE = "akima")
if(miss)
ans$z[ans$extrap] <- NA
## Error code 10 from sdsf3p indicates error code -2 from trmesh:
## first three points collinear.
## Try to add jitter to data locations to avoid collinearities,
## start with diff(range(x))*jitter*jitter.trials^1.5 and repeat for
## jitter.trials steps until success (ier=0)
if(ans$ier==10){
warning("collinear points, trying to add some jitter to avoid colinearities!")
jitter.trials <- 1
success <- FALSE
while(jitter.trials<jitter.iter & !success){
## dont use random jitter for reproducabilty by default:
if(jitter.random){
xj <- x+rnorm(n,0,diff(range(x))*jitter*jitter.trials^1.5)
yj <- y+rnorm(n,0,diff(range(y))*jitter*jitter.trials^1.5)
} else {
xj <- x+rep(c(-1,0,1),length.out=length(x))*diff(range(x))*jitter*jitter.trials^1.5
yj <- y+rep(c(0,1,-1),length.out=length(y))*diff(range(y))*jitter*jitter.trials^1.5
}
ans <- .Fortran("sdsf3p",
as.integer(1),
as.integer(n),
as.double(xj),
as.double(yj),
as.double(z),
as.integer(nx),
x = as.double(xo),
as.integer(ny),
y = as.double(yo),
z = zo,
ier = integer(1),
double(36 * n),
integer(25 * n),
extrap = as.logical(matrix(extrap,nx,ny)),
near = integer(n),
nxt = integer(n),
dist = double(n),
linear = as.logical(linear),
PACKAGE = "akima")
if(miss)
ans$z[ans$extrap] <- NA
success <- (ans$ier==0)
if(success)
warning("success: collinearities reduced through jitter")
jitter.trials <- jitter.trials+1
}
}
## prepare return value
if(is.sp){
zm <- dim(ret$z)[1]
zn <- dim(ret$z)[2]
zvec <- c(ret$z)
xvec <- c(matrix(rep(ret$x,zn),nrow=zm,ncol=zn,byrow=FALSE))
yvec <- c(matrix(rep(ret$y,zm),nrow=zm,ncol=zn,byrow=TRUE))
zm <- dim(ans$z)[1]
zn <- dim(ans$z)[2]
zvec <- c(ans$z)
xvec <- c(matrix(rep(ans$x,zn),nrow=zm,ncol=zn,byrow=FALSE))
yvec <- c(matrix(rep(ans$y,zm),nrow=zm,ncol=zn,byrow=TRUE))
nona <- !is.na(zvec)
ret <- data.frame(xvec[nona],yvec[nona],zvec[nona])
names(ret) <- c(sp.coord[1],sp.coord[2],sp.z)
coordinates(ret) <- sp.coord
ret@proj4string <- sp.proj4string
gridded(ret) <- TRUE
} else {
ret <- list(x=ans$x,y=ans$y,z=ans$z)
}
ret
}
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)
warning("'ncp' not supported, it is automatically choosen by Fortran code\n")
else
linear <- TRUE
}
if(linear)
stop("linear interpolation not implemented in interp.new().\n",
"use 'interp()' (or 'interp.old()').")
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, 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),
as.double(x),
as.double(y),
as.double(z),
as.integer(nx),
x = as.double(xo),
as.integer(ny),
y = as.double(yo),
z = zo,
ier = integer(1),
double(36 * n),
integer(25 * n),
extrap = as.logical(misso),
near = integer(n),
nxt = integer(n),
dist = double(n),
PACKAGE = "akima")[c("x", "y", "z", "extrap")]
if(miss)
ans$z[ans$extrap] <- NA
ans[c("x", "y", "z")]
}
"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)
{
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")]
}
"interpp"<-function(x, y=NULL, z, xo, yo=NULL, linear = TRUE, extrap = FALSE,
duplicate = "error", dupfun = NULL, ncp=NULL)
duplicate = "error", dupfun = NULL,
jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE)
{
# for backward compatibility
if(!is.null(ncp)){
warning('use of \'ncp\' parameter is deprecated!')
if(ncp==0)
linear <- TRUE
else if(ncp>0)
linear <- FALSE
## handle sp data, save coordinate and value names
is.sp <- FALSE
sp.coord <- NULL
sp.z <- NULL
sp.proj4string <- NULL
if(is.null(y)&&is.character(z)){
if(class(xo)=="SpatialPointsDataFrame"){
yo <- coordinates(xo)[,2]
xo <- coordinates(xo)[,1]
} else
stop(paste("either x,y,z,xo,yo have to be numeric vectors",
"or both x and xo have to be SpatialPointsDataFrames",
"and z a name of a data column in x"))
if(class(x)=="SpatialPointsDataFrame"){
sp.coord <- dimnames(coordinates(x))[[2]]
sp.z <- z
sp.proj4string <- x@proj4string
z <- x@data[,z]
y <- coordinates(x)[,2]
x <- coordinates(x)[,1]
is.sp <- TRUE
} else
stop(paste("either x,y,z,xo,yo have to be numeric vectors",
"or both x and xo have to be SpatialPointsDataFrames",
"and z a name of a data column in x"))
}
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")
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
stop('ncp < 0 ?')
}
## handle sp data, save coordinate and value names
is.sp <- FALSE
sp.coord <- NULL
sp.z <- NULL
sp.proj4string <- NULL
if(is.null(y)&&is.character(z)){
if(class(xo)