diff --git a/DESCRIPTION b/DESCRIPTION index ad2a760f11db5e8bbf66bf974dd63c9279d27d1b..7b8db178a7e883eb0921c273c5476d870b99c1ad 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,17 @@ Package: akima -Version: 0.5-11 -Date: 2013-01-19 -Title: Interpolation of irregularly spaced data +Version: 0.5-12 +Date: 2015-08-26 +Title: Interpolation of Irregularly and Regularly Spaced Data Author: Fortran code by H. Akima - R port by Albrecht Gebhardt + R port by Albrecht Gebhardt aspline function by Thomas Petzoldt interp2xyz, enhancements and corrections by Martin Maechler Maintainer: Albrecht Gebhardt -Description: Linear or cubic spline interpolation for irregular gridded data +Description: Several cubic spline interpolation methods of H. Akima for irregular and + regular gridded data are available through this package, both for the bivariate case + (irregular data: ACM 526 and ACM 761, regular data: ACM 433) and univariate case (ACM 697). + Linear interpolation of irregular gridded data is also covered by reusing D. J. Renkas + triangulation code which is part of Akimas Fortran code. License: ACM | file LICENSE Depends: R (>= 2.0.0) +Imports: sp diff --git a/NAMESPACE b/NAMESPACE index cf37c735c4cd99859a46692c3f57791b566704a4..2494840d3e5b6cd0d4a043f6cc03d2193a371dbf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,2 +1,11 @@ useDynLib(akima) + export(aspline,interp,interpp,interp.old,interp.new,interp2xyz,bicubic,bicubic.grid) + +importFrom("sp","coordinates") +importFrom("sp","coordinates<-") +importFrom("sp","gridded<-") + +importFrom("grDevices", "xy.coords") +importFrom("graphics", "hist") +importFrom("stats", "median") diff --git a/PORTING b/PORTING deleted file mode 100644 index c06b2bc339d20d88a9d9914aa755259e5dc73594..0000000000000000000000000000000000000000 --- a/PORTING +++ /dev/null @@ -1,63 +0,0 @@ -This library contains an R implementation of the S-Plus function interp(). - -See ChangeLog for details about further versions - -Version 0.1-1 - -I used the S-Plus version of interp() (S-Plus 3.3 f. Win 3.11) as starting -point. Then I searched for Akimas idsfft function, which is referenced in -interp(). - -I found the appropriate Fortran code (1978) in the ACM Collected -Algorithms archive under - - http://www.netlib.org/toms/526 - -(also included here) and splitted it into the Fortran files under src/ . -The test driver ttidbs is also included and can be compiled with "make test". - -However, it seems that this code differs a little bit from that used in S-Plus: -It implements "only" bicubic spline interpolation, no linear interpolation as -in S-Plus. - -So I modified IDSFFT and added a subroutine IDPTLI which does linear -interpolation within the triangles, generated by IDTANG. - -Further changes are -+ REAL -> DOUBLE PRECISION -+ static DIMENSIONs replaced with dynamic -+ option to toggle extrapolation outside of convex hull added in IDSFFT and - IDBVIP. - Because I don't know how to generate NAs in Fortran (I use g77 on Linux), - I added a logical array MISSI, that indicates if a returned value should - be NA. These values will be set to NA after the Fortran call. -+ option to handle duplicate data points added (according to an example in the - S-Plus help page) -+ man pages converted and rewritten -+ data set akima (from S-Plus) added -+ function interpp() added, it evaluates the interpolated function at - arbitraryly choosen points and generates no regular grid as interp() does. - There where some problems with interpp() when using the Fortran version: - - it crashes when compiled with "g77 -O2 -fpic" and called with more than - one output point. - - compilation with "g77 -g -fpic" fails (see src/Makefile for details) - - compilation with "g77 -g" works (and no crashes occur) - These problems do not occur in the C Version (generated by f2c), so it would - be better to use only the src-c tree. - - - -After I finished the above steps I found a more recent version (ACM 761, 1996) -of Akimas interpolation code which uses the tripack package (also available at -ACM as algorithm no. 751) for triangulation and now I'm trying to use it for -the next version of interp(). - - ------------------------------------------------------------------- -Albrecht Gebhardt email: albrecht.gebhardt@uni-klu.ac.at -Institut fuer Mathematik Tel. : (++43 463) 2700/837 -Universitaet Klagenfurt Fax : (++43 463) 2700/834 -Villacher Str. 161 -A-9020 Klagenfurt, Austria ------------------------------------------------------------------- - diff --git a/R/interp.R b/R/interp.R index 7838733621820bf82316bdf4934e7b528f4ea32a..39eec3571d7f9965f1e62be7d68eaea4dfee2de4 100644 --- a/R/interp.R +++ b/R/interp.R @@ -1,26 +1,64 @@ interp <- - function(x, y, z, - xo = seq(min(x), max(x), length = 40), - yo = seq(min(y), max(y), length = 40), linear=TRUE, - extrap = FALSE, duplicate = "error", dupfun = NULL, ncp=NULL) + 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) { - # 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 ?') - } - if(linear) - ## use the old version for linear interpolation - interp.old(x, y, z, xo = xo, yo = yo, ncp = 0, - 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) + ## 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 + sp.z <- NULL + sp.proj4string <- NULL + if(is.null(y)&&is.character(z)){ + 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 + xo = seq(min(x), max(x), length = nx) + yo = seq(min(y), max(y), length = ny) + } else + 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(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)) + 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 + } + ret } diff --git a/R/interpp.R b/R/interpp.R index ad2d202fc17060eb57a27c7c564b93f4b5f1ce93..7cfb8ab747dd4252ea82cbeced361d45a9b06b1e 100644 --- a/R/interpp.R +++ b/R/interpp.R @@ -1,4 +1,4 @@ -"interpp"<-function(x, y, z, xo, yo, linear = TRUE, extrap = FALSE, +"interpp"<-function(x, y=NULL, z, xo, yo=NULL, linear = TRUE, extrap = FALSE, duplicate = "error", dupfun = NULL, ncp=NULL) { @@ -12,11 +12,43 @@ 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)=="SpatialPointsDataFrame"){ + yo <- coordinates(xo)[,2] + xo <- coordinates(xo)[,1] + } else + stop("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("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(linear) - # the old Akima code: - interpp.old(x, y, z, xo, yo, ncp=0, extrap, duplicate, dupfun) + ## the old Akima code: + ret <- interpp.old(x, y, z, xo, yo, ncp=0, extrap, duplicate, dupfun) else - # new code for splines - interpp.new(x, y, z, xo, yo, extrap, duplicate, dupfun) + ## new code for splines + ret <- interpp.new(x, y, z, xo, yo, extrap, duplicate, dupfun) + if(is.sp){ + nona <- !is.na(ret$z) + ret <- data.frame(ret$x[nona],ret$y[nona],ret$z[nona]) + names(ret) <- c(sp.coord[1],sp.coord[2],sp.z) + coordinates(ret) <- sp.coord + ret@proj4string <- sp.proj4string + } + ret } diff --git a/README b/README index 3fbfc9d564be3ae44677d15045c043eaee9d8784..44c555eedcfc72b940830af2e9c0765037e0a1c1 100644 --- a/README +++ b/README @@ -21,11 +21,83 @@ Thomas Petzoldt provided an R interface to Akimas univariate spline interpolation functions (aspline) which has been included into this library. +The original papers from H. Akima are available (access is not free, depends on your +ACM subscription): +TOMS 433: http://dl.acm.org/citation.cfm?id=355604.355605 +TOMS 526: http://dl.acm.org/citation.cfm?id=355780.355787 +TOMS 697: http://dl.acm.org/citation.cfm?id=114697.116808 +TOMS 761: http://dl.acm.org/citation.cfm?id=232826.232856 + +Original Fortran code is freely available at Netlib/TOMS -- Collected Algorithms +of the ACM at: +http://www.netlib.org/toms/433 +http://www.netlib.org/toms/526 +http://www.netlib.org/toms/697 +http://www.netlib.org/toms/761 + ------------------------------------------------------------------------------ -Albrecht Gebhardt email: albrecht.gebhardt@uni-klu.ac.at +Albrecht Gebhardt email: albrecht.gebhardt@aau.at Institut fuer Mathematik Tel. : (++43 463) 2700/3118 Universitaet Klagenfurt Fax : (++43 463) 2700/3198 -Universitaetsstr. 65-67 WWW : http://www.math.uni-klu.ac.at/~agebhard +Universitaetsstr. 65-67 WWW : http://www.stat.aau.at/~agebhard A-9020 Klagenfurt, Austria ------------------------------------------------------------------------------ + + + +Details: + +See ChangeLog for details about further versions + +Version 0.1-1 + +I used the S-Plus version of interp() (S-Plus 3.3 f. Win 3.11) as starting +point. Then I searched for Akimas idsfft function, which is referenced in +interp(). + +I found the appropriate Fortran code (1978) in the ACM Collected +Algorithms archive under + + http://www.netlib.org/toms/526 + +(also included here) and splitted it into the Fortran files under src/ . +The test driver ttidbs was also included in older versions and could be +compiled with "make test", now it is removed here. + +However, it seems that this code differs a little bit from that used in S-Plus: +It implements "only" bicubic spline interpolation, no linear interpolation as +in S-Plus. + +So I modified IDSFFT and added a subroutine IDPTLI which does linear +interpolation within the triangles, generated by IDTANG. + +Further changes are ++ REAL -> DOUBLE PRECISION ++ static DIMENSIONs replaced with dynamic ++ option to toggle extrapolation outside of convex hull added in IDSFFT and + IDBVIP. + Because I don't know how to generate NAs in Fortran (I use g77 on Linux), + I added a logical array MISSI, that indicates if a returned value should + be NA. These values will be set to NA after the Fortran call. ++ option to handle duplicate data points added (according to an example in the + S-Plus help page) ++ man pages converted and rewritten ++ data set akima (from S-Plus) added ++ function interpp() added, it evaluates the interpolated function at + arbitraryly choosen points and generates no regular grid as interp() does. + There where some problems with interpp() when using the Fortran version: + - it crashes when compiled with "g77 -O2 -fpic" and called with more than + one output point. + - compilation with "g77 -g -fpic" fails (see src/Makefile for details) + - compilation with "g77 -g" works (and no crashes occur) + These problems do not occur in the C Version (generated by f2c), so it would + be better to use only the src-c tree. + + + +After I finished the above steps I found a more recent version (ACM 761, 1996) +of Akimas interpolation code which uses the tripack package (also available at +ACM as algorithm no. 751) for triangulation and now I'm trying to use it for +the next version of interp(). + diff --git a/TITLE b/TITLE deleted file mode 100644 index bf92d1d4381b606b7d53b2d6425340ee1660f562..0000000000000000000000000000000000000000 --- a/TITLE +++ /dev/null @@ -1 +0,0 @@ -akima Interpolation of irregularly spaced data diff --git a/cleanup b/cleanup index 31686c7de75af08ab6242fd59ee20603a8c4dbc9..21a17384193fdbaa7ed0ea63cf767bdafa1e4ba0 100755 --- a/cleanup +++ b/cleanup @@ -4,4 +4,5 @@ rm -f src/Makevars rm -f src/*.o rm -f src/*.so rm -f src/*.dll +rm -f config.* diff --git a/man/interp.Rd b/man/interp.Rd index eae204e1336c23dce2b165b3db1e91d76b5e4600..cad82a5f7952d5822202f34894ba3b7d3ac694bc 100644 --- a/man/interp.Rd +++ b/man/interp.Rd @@ -10,9 +10,10 @@ Akima. } \usage{ -interp(x, y, z, xo=seq(min(x), max(x), length = 40), - yo=seq(min(y), max(y), length = 40), - linear = TRUE, extrap=FALSE, duplicate = "error", dupfun = NULL, ncp = NULL) +interp(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) 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) @@ -22,18 +23,27 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40), } \arguments{ \item{x}{ - vector of x-coordinates of data points. + 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. + 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 and may + \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 @@ -59,7 +69,7 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40), deprecated, use parameter \code{linear}. Now only used by \code{interp.old()}. - meaning was: + 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 @@ -80,6 +90,8 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40), }} \item{dupfun}{a function, applied to duplicate points if \code{duplicate= "user"}.} + \item{nx}{dimension of output grid in x direction} + \item{ny}{dimension of output grid in y direction} } \value{ list with 3 components: @@ -92,6 +104,9 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40), 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} is a wrapper for the two versions \code{interp.old} (it @@ -163,12 +178,12 @@ akima.smooth <- image (akima.smooth, main = "interp(, *) on finer grid") contour(akima.smooth, add = TRUE, col = "thistle") points(akima, pch = 3, cex = 2, col = "blue") -# use triangulation package to show underlying triangulation: +## use triangulation package to show underlying triangulation: \dontrun{ if(library(tripack, logical.return=TRUE)) plot(tri.mesh(akima), add=TRUE, lty="dashed") } -# use only 15 points (interpolation only within convex hull!) +## use only 15 points (interpolation only within convex hull!) akima.part <- with(akima, interp(x[1:15], y[1:15], z[1:15])) image(akima.part) title("interp() on subset of only 15 points") @@ -206,20 +221,32 @@ filled.contour(akima.spl, color.palette = full.pal, plot.axes = { axis(1); axis(2); title("smooth interp(*, linear = FALSE)"); points(akima, pch = 3, col= hcl(c=100, l = 20))}) -# no extrapolation! +## no extrapolation! ## example with duplicate points : data(airquality) air <- subset(airquality, !is.na(Temp) & !is.na(Ozone) & !is.na(Solar.R)) -# gives an error {duplicate ..}: +## gives an error {duplicate ..}: try( air.ip <- interp(air$Temp,air$Solar.R,air$Ozone, linear=FALSE) ) -# use mean of duplicate points: +## use mean of duplicate points: air.ip <- with(air, interp(Temp, Solar.R, log(Ozone), duplicate = "mean", linear = FALSE)) image(air.ip, main = "Airquality: Ozone vs. Temp and Solar.R") with(air, points(Temp, Solar.R)) + +\dontrun{ + ## interp can handle spatial point dataframes created by the sp package: + library(sp) + data(meuse) + coordinates(meuse) <- ~x+y + ## argument z has to be named, y has to be omitted! + z <- interp(meuse,z="zinc",nx=100,ny=150) + spplot(z,"zinc") + z <- interp(meuse,z="zinc",nx=100,ny=150,linear=FALSE) + spplot(z,"zinc") +} } \keyword{dplot} diff --git a/man/interpp.Rd b/man/interpp.Rd index 0b0ffbc4a7a36eaa484d7d4787c0bbf7a61e1024..278b336737f6d7e3ebe7db35fc9a5f3cd0f4a014 100644 --- a/man/interpp.Rd +++ b/man/interpp.Rd @@ -6,43 +6,57 @@ \alias{interpp.old} \alias{interpp.new} \usage{ -interpp(x, y, z, xo, yo, linear=TRUE, extrap=FALSE, duplicate = "error", -dupfun = NULL, ncp) +interpp(x, y=NULL, z, xo, yo=NULL, linear=TRUE, extrap=FALSE, +duplicate = "error", dupfun = NULL, ncp) } \arguments{ \item{x}{ - vector of x-coordinates of data points. + 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. + 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 and may contain no fewer + + \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). + accepted). } \item{xo}{ vector of x-coordinates of points at which to evaluate the interpolating - function.} + 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.} + function. + + If operating on \code{SpatialPointsDataFrame}s this is left as \code{NULL} + } \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 + 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. + 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). @@ -57,12 +71,12 @@ dupfun = NULL, ncp) calculate mean , median or user defined function of duplicate z values. } -\item{dupfun}{this function is applied to duplicate points if \code{duplicate="user"} -} +\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}. @@ -74,35 +88,37 @@ dupfun = NULL, ncp) \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} refer to Akimas Fortran code from 1978 and 1996 resp. The call wrapper \code{interpp} - chooses \code{interpp.old} for linear and \code{interpp.new} for cubic + chooses \code{interpp.old} for linear and \code{interpp.new} for cubic spline interpolation. - Earlier versions (pre 0.5-1) of \code{interpp} used the parameter + 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 + 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. + 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}. + 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 @@ -120,8 +136,8 @@ dupfun = NULL, ncp) \bold{22}, 362-371. } \seealso{ - \code{\link[graphics]{contour}}, \code{\link[graphics]{image}}, - \code{\link[stats]{approxfun}}, \code{\link[stats]{splinefun}}, + \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}}. } @@ -134,5 +150,19 @@ akima.lip$z akima.sip<-interpp(akima$x, akima$y, akima$z,c(1,5,10),c(2,6,12), linear=FALSE) akima.sip$z +\dontrun{ + ## interaction with sp objects: + library(sp) + ## take 30 sample points out of meuse grid: + data(meuse.grid) + m0 <- meuse.grid[sample(1:3103,30),] + coordinates(m0) <- ~x+y + ## interpolate on this 30 points: + ## note: both "meuse" and "m0" are sp objects + ## (SpatialPointsDataFrame) !! + ## arguments z and xo have to named, y has to be omitted! + ipp <- interpp(meuse,z="zinc",xo=m0) + spplot(ipp) +} } \keyword{dplot} diff --git a/src/ttidbs.f b/src/ttidbs.f deleted file mode 100644 index 7d594d8544446731b90359e456f9b72e010c126f..0000000000000000000000000000000000000000 --- a/src/ttidbs.f +++ /dev/null @@ -1,124 +0,0 @@ - PROGRAM TTIDBS -C PROGRAM TTIDBS(OUTPUT,TAPE6=OUTPUT) ID000070 -C THIS PROGRAM IS A TEST PROGRAM FOR THE IDBVIP/IDSFFT SUBPRO- ID000080 -C GRAM PACKAGE. ALL ELEMENTS OF RESULTING DZI1 AND DZI2 ARRAYS ID000090 -C ARE EXPECTED TO BE ZERO. ID000100 -C THE LUN CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS ID000110 -C THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, ID000120 -C THEREFORE, SYSTEM DEPENDENT. ID000130 -C DECLARATION STATEMENTS ID000140 - IMPLICIT DOUBLE PRECISION (A-D,P-Z) - DIMENSION XD(30),YD(30),ZD(30), ID000150 - 1 XI(6),YI(5),ZI(6,5),MISSI(6,5), ID000160 - 2 ZI1(6,5),ZI2(6,5),DZI1(6,5),DZI2(6,5), ID000170 - 3 IWK(1030),WK(240) ID000180 - LOGICAL MISSI - DATA NCP/4/ ID000190 - DATA NDP/30/ ID000200 - DATA XD(1), XD(2), XD(3), XD(4), XD(5), XD(6), ID000210 - 1 XD(7), XD(8), XD(9), XD(10),XD(11),XD(12), ID000220 - 2 XD(13),XD(14),XD(15),XD(16),XD(17),XD(18), ID000230 - 3 XD(19),XD(20),XD(21),XD(22),XD(23),XD(24), ID000240 - 4 XD(25),XD(26),XD(27),XD(28),XD(29),XD(30)/ ID000250 - 5 11.16, 24.20, 19.85, 10.35, 19.72, 0.00, ID000260 - 6 20.87, 19.99, 10.28, 4.51, 0.00, 16.70, ID000270 - 7 6.08, 25.00, 14.90, 0.00, 9.66, 5.22, ID000280 - 8 11.77, 15.10, 25.00, 25.00, 14.59, 15.20, ID000290 - 9 5.23, 2.14, 0.51, 25.00, 21.67, 3.31/ ID000300 - DATA YD(1), YD(2), YD(3), YD(4), YD(5), YD(6), ID000310 - 1 YD(7), YD(8), YD(9), YD(10),YD(11),YD(12), ID000320 - 2 YD(13),YD(14),YD(15),YD(16),YD(17),YD(18), ID000330 - 3 YD(19),YD(20),YD(21),YD(22),YD(23),YD(24), ID000340 - 4 YD(25),YD(26),YD(27),YD(28),YD(29),YD(30)/ ID000350 - 5 1.24, 16.23, 10.72, 4.11, 1.39, 20.00, ID000360 - 6 20.00, 4.62, 15.16, 20.00, 4.48, 19.65, ID000370 - 7 4.58, 11.87, 3.12, 0.00, 20.00, 14.66, ID000380 - 8 10.47, 17.19, 3.87, 0.00, 8.71, 0.00, ID000390 - 9 10.72, 15.03, 8.37, 20.00, 14.36, 0.13/ ID000400 - DATA ZD(1), ZD(2), ZD(3), ZD(4), ZD(5), ZD(6), ID000410 - 1 ZD(7), ZD(8), ZD(9), ZD(10),ZD(11),ZD(12), ID000420 - 2 ZD(13),ZD(14),ZD(15),ZD(16),ZD(17),ZD(18), ID000430 - 3 ZD(19),ZD(20),ZD(21),ZD(22),ZD(23),ZD(24), ID000440 - 4 ZD(25),ZD(26),ZD(27),ZD(28),ZD(29),ZD(30)/ ID000450 - 5 22.15, 2.83, 7.97, 22.33, 16.83, 34.60, ID000460 - 6 5.74, 14.72, 21.59, 15.61, 61.77, 6.31, ID000470 - 7 35.74, 4.40, 21.70, 58.20, 4.73, 40.36, ID000480 - 8 13.62, 12.57, 8.74, 12.00, 14.81, 21.60, ID000490 - 9 26.50, 53.10, 49.43, 0.60, 5.52, 44.08/ ID000500 - DATA NXI/6/, NYI/5/ ID000510 - DATA XI(1), XI(2), XI(3), XI(4), XI(5), XI(6)/ ID000520 - 1 0.00, 5.00, 10.00, 15.00, 20.00, 25.00/ ID000530 - DATA YI(1), YI(2), YI(3), YI(4), YI(5)/ ID000540 - 1 0.00, 5.00, 10.00, 15.00, 20.00/ ID000550 - DATA ZI(1,1),ZI(2,1),ZI(3,1),ZI(4,1),ZI(5,1),ZI(6,1), ID000560 - 1 ZI(1,2),ZI(2,2),ZI(3,2),ZI(4,2),ZI(5,2),ZI(6,2), ID000570 - 2 ZI(1,3),ZI(2,3),ZI(3,3),ZI(4,3),ZI(5,3),ZI(6,3), ID000580 - 3 ZI(1,4),ZI(2,4),ZI(3,4),ZI(4,4),ZI(5,4),ZI(6,4), ID000590 - 4 ZI(1,5),ZI(2,5),ZI(3,5),ZI(4,5),ZI(5,5),ZI(6,5)/ ID000600 - 5 58.20, 39.55, 26.90, 21.71, 17.68, 12.00, ID000610 - 6 61.58, 39.39, 22.04, 21.29, 14.36, 8.04, ID000620 - 7 59.18, 27.39, 16.78, 13.25, 8.59, 5.36, ID000630 - 8 52.82, 40.27, 22.76, 16.61, 7.40, 2.88, ID000640 - 9 34.60, 14.05, 4.12, 3.17, 6.31, 0.60/ ID000650 - DATA LUN/6/ ID000660 -C CALCULATION ID000670 - 10 MD=1 ID000680 - DO 12 IYI=1,NYI ID000690 - DO 11 IXI=1,NXI ID000700 - IF(IXI.NE.1.OR.IYI.NE.1) MD=2 ID000710 - MISSI(IXI,IYI)=.FALSE. - CALL IDBVIP(MD,NCP,NDP,XD,YD,ZD,1,XI(IXI),YI(IYI), ID000720 - 1 ZI1(IXI,IYI),IWK,WK,MISSI) ID000730 - 11 CONTINUE ID000740 - 12 CONTINUE ID000750 - 15 CALL IDSFFT(1,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI2,IWK,WK,MISSI) ID000760 - DO 17 IYI=1,NYI ID000770 - DO 16 IXI=1,NXI ID000780 - DZI1(IXI,IYI)=ABS(ZI1(IXI,IYI)-ZI(IXI,IYI)) ID000790 - DZI2(IXI,IYI)=ABS(ZI2(IXI,IYI)-ZI(IXI,IYI)) ID000800 - 16 CONTINUE ID000810 - 17 CONTINUE ID000820 -C PRINTING OF INPUT DATA ID000830 -C WRITE (LUN,2020) NDP ID000840 - DO 23 IDP=1,NDP ID000850 -C WRITE (LUN,2021) ID000860 -C WRITE (LUN,2022) IDP,XD(IDP),YD(IDP),ZD(IDP) ID000870 - 23 CONTINUE ID000880 -C PRINTING OF OUTPUT RESULTS ID000890 -C WRITE (LUN,2030) ID000900 -C WRITE (LUN,2031) YI ID000910 - DO 33 IXI=1,NXI ID000920 -C WRITE (LUN,2032) XI(IXI),(ZI1(IXI,IYI),IYI=1,NYI) ID000930 - 33 CONTINUE ID000940 -C WRITE (LUN,2040) ID000950 -C WRITE (LUN,2031) YI ID000960 - DO 43 IXI=1,NXI ID000970 -C WRITE (LUN,2032) XI(IXI),(DZI1(IXI,IYI),IYI=1,NYI) ID000980 - 43 CONTINUE ID000990 -C WRITE (LUN,2050) ID001000 -C WRITE (LUN,2031) YI ID001010 - DO 53 IXI=1,NXI ID001020 -C WRITE (LUN,2032) XI(IXI),(ZI2(IXI,IYI),IYI=1,NYI) ID001030 - 53 CONTINUE ID001040 -C WRITE (LUN,2060) ID001050 -C WRITE (LUN,2031) YI ID001060 - DO 63 IXI=1,NXI ID001070 -C WRITE (LUN,2032) XI(IXI),(DZI2(IXI,IYI),IYI=1,NYI) ID001080 - 63 CONTINUE ID001090 -C STOP ID001100 -C FORMAT STATEMENTS ID001110 - 2020 FORMAT(1H1,6HTTIDBS/////3X,10HINPUT DATA,8X,5HNDP =,I3/// ID001120 - 1 30H I XD YD ZD /) ID001130 - 2021 FORMAT(1X) ID001140 - 2022 FORMAT(5X,I2,2X,3F7.2) ID001150 - 2030 FORMAT(1H1,6HTTIDBS/////3X,17HIDBVIP SUBROUTINE/// ID001160 - 1 26X,10HZI1(XI,YI)) ID001170 - 2031 FORMAT(7X,2HXI,4X,3HYI=/12X,5F7.2/) ID001180 - 2032 FORMAT(1X/1X,F9.2,2X,5F7.2) ID001190 - 2040 FORMAT(1X/////3X,10HDIFFERENCE/// ID001200 - 1 25X,11HDZI1(XI,YI)) ID001210 - 2050 FORMAT(1H1,6HTTIDBS/////3X,17HIDSFFT SUBROUTINE/// ID001220 - 1 26X,10HZI2(XI,YI)) ID001230 - 2060 FORMAT(1X/////3X,10HDIFFERENCE/// ID001240 - 1 25X,11HDZI2(XI,YI)) ID001250 - END ID001260