From 13677bba35b0aa073f967964f4a797abad4f1b0e Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 17 Nov 2017 12:46:59 +0100 Subject: [PATCH 01/14] remove Trackstat.R to test --- R/Trackstat.R | 412 -------------------------------------------------- 1 file changed, 412 deletions(-) delete mode 100644 R/Trackstat.R diff --git a/R/Trackstat.R b/R/Trackstat.R deleted file mode 100644 index 7f1010d..0000000 --- a/R/Trackstat.R +++ /dev/null @@ -1,412 +0,0 @@ -# function as.Track accepts a dataframe (let say X) with three columns "xcoor", "ycoor" and time (or any class that can be converted to a data.frame) -# and converts it to an object of class Track. It can also accepts covariates for the corresponding -# locations, covariates must be a dataframe with some columns and length of each column is equal -# to length of each column in X. -as.Track <- function(X,covariate){ - stopifnot(nrow(X)>0) - colnames(X) <- c("xcoor","ycoor","time") - if(!is.data.frame(X)) X <- as.data.frame(X) - sp <- cbind(x=X$xcoor,y=X$ycoor) - sp <- SpatialPoints(sp) - t <- as.POSIXct(paste(X$date,X$time)) - if(missing(covariate)) covariate <- data.frame(d=rep(NA,length(X$xcoor))) - - return(Track(STIDF(sp,time = t,data =covariate))) -} - - -# function reTrack accepts X as an object of class Track. Output is a reconstructed Track (an object of class Track), based on "timestamp". -# It only returns the interpolated points. -reTrack <- function(X,at=c("track","dfrm"),timestamp=timestamp,tsq=NULL){ - - if (missing(tsq)) tsq <- tsqTracks(X,timestamp = timestamp) - if(missing(at)) at <- "track" - Xrange <- rngTrack(X) - X <- cbind(as.data.frame(X)[c(coordnames(X), "time")]) - xnew <- c() - ynew <- c() - - x <- X$x - y <- X$y - - time <- tsq[tsqXrange[1]] - ivs <- findInterval(time,X$time) - - for (i in 1:length(ivs)) { - if (!ivs[i] == 0 && !ivs[i] == nrow(X)) { - - iv <- ivs[i] - tdiff1 <- difftime(time[i],X$time[iv],units = "sec") # diff between timestamp and start of the interval it falls in - tdiff2 <- difftime(X$time[iv+1],X$time[iv],units = units(tdiff1)) # diff between timestamps (calculated here because it often varies) - ratio <- as.numeric(tdiff1)/as.numeric(tdiff2) - x1 <- X[iv,1] # segment coordinates - y1 <- X[iv,2] - x2 <- X[iv+1,1] - y2 <- X[iv+1,2] - xnew <- c(xnew, x1 + ratio * (x2 - x1)) #find point - ynew <- c(ynew, y1 + ratio * (y2 - y1)) - } - } - - newTrack <- data.frame(xnew, ynew, time) - newTrack <- newTrack[!duplicated(newTrack),] # remove duplicates - newTrack <- newTrack[order(newTrack$time),] # sort by timestamp - colnames(newTrack) <- c("xcoor","ycoor","time") - if (at=="dfrm") {attr(newTrack,"tsq") <-tsq;return(newTrack) } - return(as.Track(newTrack)) -} - -# rngTrack returns the timerange of an object of class Track -rngTrack <- function(X){ - Y <- cbind(as.data.frame(X)[c(coordnames(X), "time")]) - return(range(Y$time)) -} - -# tsqtracks returns a sequance of time based on a list of tracks (or a single object of class Track) and an argument timestamp -tsqTracks <- function(X,timestamp){ - - if (!is.list(X)) timerange <- rngTrack(X) - else timerange <- lapply(X, rngTrack) - - Trackrg <- range(timerange) - class(Trackrg) <- c('POSIXt','POSIXct') - # a seq from the range has been created every timestamp - timeseq <- seq(from=as.POSIXct(strftime(Trackrg[1])),to=as.POSIXct(strftime(Trackrg[2])),by = timestamp) - - return(timeseq) - -} - - - -# function avedistTrack accepts X as a list of tracks and reports the average distance between -# tracks over time, output is an object of class "distrack" -avedistTrack <- function(X,timestamp){ - - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - # calculate a sequance of time to interpolate tracks within this sequance - timeseq <- tsqTracks(X,timestamp = timestamp) - - # reconstruct tracks in sequance timeseq - Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") - - wincor <- lapply(X=1:length(Z),FUN = function(i){ - return(list(min(Z[[i]]$x),max(Z[[i]]$x),min(Z[[i]]$y),max(Z[[i]]$y))) - }) - wincor <- matrix(unlist(wincor),nrow = 4) - w <- owin(c(min(wincor[1,]),max(wincor[2,])),c(min(wincor[3,]),max(wincor[4,]))) - # create a list and convert tracks in each element of timeseq to an object of calss ppp - p <- list() - - for (j in 1:length(timeseq)) { - - l <- lapply(X=1:length(Z), function(i){ - Z[[i]][which(Z[[i]]$time==timeseq[j]),-3] - }) - - if(length(unlist(l))>0){ - x <- unlist(lapply(X=1:length(l),function(k){ - l[[k]]$x - })) - y <- unlist(lapply(X=1:length(l),function(h){ - l[[h]]$y - })) - p[[j]] <- as.ppp(data.frame(x,y),W=w) - } - - } - - p1 <- p[!sapply(p, is.null)] - avedist <- unlist( - lapply(X=1:length(p1), function(i){ - pd <- pairdist(p1[[i]]) - pd <- pd[pd>0] - return(mean(pd)) - }) - ) - avedist <- data.frame(timeseq[!sapply(p, is.null)],avedist) - class(avedist) <- c("distrack") - attr(avedist,"ppp") <- p - return(avedist) -} - -print.distrack <- function(x){ - print(as.vector(x$avedist)) -} - -plot.distrack <- function(x,...){ - plot(x$timeseq,x$avedist,,xlab="time",ylab="average distance",...) -} - - -rmdupTrack <- function(X){ - X <- cbind(as.data.frame(X)[c(coordnames(X), "time")]) - X <- X[!duplicated(X),] - return(as.Track(X)) -} - -as.Track.ppp <- function(X,timestamp){ - - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - # calculate a sequance of time to interpolate tracks within this sequance - timeseq <- tsqTracks(X,timestamp = timestamp) - - # reconstruct tracks in sequance timeseq - Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") - - wincor <- lapply(X=1:length(Z),FUN = function(i){ - return(list(min(Z[[i]]$x),max(Z[[i]]$x),min(Z[[i]]$y),max(Z[[i]]$y))) - }) - wincor <- matrix(unlist(wincor),nrow = 4) - w <- owin(c(min(wincor[1,]),max(wincor[2,])),c(min(wincor[3,]),max(wincor[4,]))) - - # create a list and convert tracks in each element of timeseq to an object of calss ppp - p <- list() - - for (j in 1:length(timeseq)) { - - l <- lapply(X=1:length(Z), function(i){ - w <- which(Z[[i]]$time==timeseq[j]) - if (length(w)>0) return(cbind(Z[[i]][w,-3],id=i)) - return(Z[[i]][w,-3]) - }) - - if(length(unlist(l))>0){ - x <- unlist(lapply(X=1:length(l),function(k){ - l[[k]]$x - })) - y <- unlist(lapply(X=1:length(l),function(h){ - l[[h]]$y - })) - m <- unlist(lapply(X=1:length(l),function(h){ - l[[h]]$id - })) - p[[j]] <- as.ppp(data.frame(x,y),W=w) - marks(p[[j]]) <- m - } - - } - return(p) -} - - -density.Track <- function(X,timestamp,...){ - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - - p <- as.Track.ppp(X,timestamp) - p <- p[!sapply(p, is.null)] - imlist <- lapply(p, density.ppp,...) - out <- Reduce("+",imlist)/length(imlist) - attr(out,"Tracksim") <- imlist - attr(out,"ppps") <- p - return(out) -} - -as.Track.arrow <- function(X,timestamp,epsilon=epsilon){ - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - - Z <- as.Track.ppp(X,timestamp) - Z <- Z[!sapply(Z, is.null)] - wind <- Z[[1]]$window - arrows <- list() - Y <- list() - for (i in 1:length(Z)) { - if(i==length(Z)) break() - j <- i+1 - m1 <- match(marks(Z[[i]]),marks(Z[[j]])) - m2 <- match(marks(Z[[j]]),marks(Z[[i]])) - m1 <- m1[!is.na(m1)] - m2 <- m2[!is.na(m2)] - x <- Z[[j]][m1] - y <- Z[[i]][m2] - l <- psp(y$x,y$y,x$x,x$y,window = wind) - arrows[[i]] <- l - center <- midpoints.psp(l) - mark <- lengths.psp(l) - marks(center) <- mark - if (missing(epsilon)) epsilon <- 0 - Y[[i]] <- center[mark>epsilon] - } - class(Y) <- c("list","Trrow") - attr(Y,"psp") <- arrows - return(Y) -} - -print.Trrow <- function(x, ...) { - attributes(x) <- NULL - print(x) -} - -Track.idw <- function(X,timestamp,epsilon=epsilon,...){ - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - - Y <- as.Track.arrow(X,timestamp,epsilon=epsilon) - Z <- lapply(Y, idw,...) - meanIDW <- Reduce("+",Z)/length(Z) - return(meanIDW) -} - -avemove <- function(X,timestamp,epsilon=epsilon){ - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - timeseq <- tsqTracks(X,timestamp = timestamp) - Y <- as.Track.arrow(X,timestamp,epsilon=epsilon) - Z <- attr(Y,"psp") - preout <- lapply(X=1:length(Z), function(i){ - mean(lengths.psp(Z[[i]])) - }) - out <- unlist(preout) - class(out) <- c("numeric", "arwlen") - attr(out,"tsq") <- timeseq - return(out) -} - -print.arwlen <- function(x){ - print(as.vector(x)) -} - -plot.arwlen <- function(x,...){ - tsq <- attr(x,"tsq") - tsq <- tsq[c(-1,-length(tsq))] - plot(tsq,x,xlab="time",ylab="average movement",...) -} - -chimaps <- function(X,timestamp,rank,...){ - if (!is.numeric(rank)) stop("rank must be numeric") - if (rank < 1 | rank >length(X)) stop("rank must be number between one and the number of Tracks") - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - timeseq <- tsqTracks(X,timestamp = timestamp) - d <- density.Track(X,timestamp,...) - imlist <- attr(d,"Tracksim") - sumim <- Reduce("+",imlist) - chi <- lapply(X=1:length(imlist),FUN = function(i){ - E1 <- sumim*sum(imlist[[i]]$v)/(sum(sumim$v)) - return((imlist[[i]]-E1)/sqrt(E1)) - }) - out <- chi[[rank]] - attr(out,"ims") <- chi - attr(out,"time") <- timeseq[rank] - attr(out,"timevec") <- timeseq - return(out) -} - -Kinhom.Track <- function(X,timestamp, - correction=c("border", "bord.modif", "isotropic", "translate"),q,...){ - - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - - cor <- match.arg(correction,correction) - - ZZ <- density.Track(X,timestamp) - - Z <- attr(ZZ,"Tracksim") - Y <- attr(ZZ,"ppps") - W <- Y[[1]]$window - ripley <- min(diff(W$xrange), diff(W$yrange))/4 - rr <- seq(0,ripley,length.out = 513) - - K <- lapply(X=1:length(Y), function(i){ - kk <- Kinhom(Y[[i]],lambda = Z[[i]],correction=cor,r=rr,...) - return(as.data.frame(kk)) - }) - Kmat <- matrix(nrow = length(K[[1]]$theo),ncol = length(K)) - for (i in 1:length(K)) { - Kmat[,i] <- K[[i]][,3] - } - # Kmat <- as.data.frame(K) - lowk <- numeric() - upk <- numeric() - avek <- numeric() - for (i in 1:nrow(Kmat)) { - avek[i] <- mean(Kmat[i,]) - lowk[i] <- quantile(Kmat[i,],q) - upk[i] <- quantile(Kmat[i,],1-q) - } - - out <- data.frame(lowk=lowk,upk=upk,avek=avek,r=K[[1]]$r,theo=K[[1]]$theo) - class(out) <- c("list","KTrack") - attr(out,"out") <- out - return(out) -} -print.KTrack <- function(x){ - print("variability area of K-function") -} - -plot.KTrack <- function(x,type="l",col= "grey70",...){ - ylim <- c(min(x$lowk),max(x$upk)) - plot(x$r,x$lowk,ylim=ylim,xlab="r",ylab="K(r)",type=type,...) - points(x$r,x$upk,type=type) - polygon(c(x$r, rev(x$r)), c(x$upk, rev(x$lowk)), - col = col, border = NA) - points(x$r,x$theo,type=type,col=2) - points(x$r,x$avek,type=type) - legend(0,max(x$upk),col = c(2,1),legend=c("poisson","average"),lty=c(1,1)) -} - -pcfinhom.Track <- function(X,timestamp, - correction = c("translate", "Ripley"),q,...){ - - stopifnot(length(X)>1 & is.list(X)) - - if (missing(timestamp)) stop("set timestamp") - - cor <- match.arg(correction,correction) - - ZZ <- density.Track(X,timestamp) - - Z <- attr(ZZ,"Tracksim") - Y <- attr(ZZ,"ppps") - - g <- lapply(X=1:length(Y), function(i){ - gg <- pcfinhom(Y[[i]],lambda = Z[[i]],correction=cor,...) - return(as.data.frame(gg)) - }) - gmat <- matrix(nrow = length(g[[1]]$theo),ncol = length(g)) - for (i in 1:length(g)) { - gmat[,i] <- g[[i]][,3] - } - # Kmat <- as.data.frame(K) - lowg <- numeric() - upg <- numeric() - aveg <- numeric() - for (i in 1:nrow(gmat)) { - aveg[i] <- mean(gmat[i,]) - lowg[i] <- quantile(gmat[i,],q) - upg[i] <- quantile(gmat[i,],1-q) - } - - out <- data.frame(lowg=lowg,upg=upg,aveg=aveg,r=g[[1]]$r,theo=g[[1]]$theo) - class(out) <- c("list","gTrack") - attr(out,"out") <- out - return(out) -} - - -print.gTrack <- function(x){ - print("variability area of pair correlatio function") -} - -plot.gTrack <- function(x,type="l",col= "grey70",...){ - ylim <- c(min(x$lowg),max(x$upg)) - plot(x$r,x$lowg,ylim=ylim,xlab="r",ylab="g",type=type,...) - points(x$r,x$upg,type=type) - polygon(c(x$r, rev(x$r)), c(x$upg, rev(x$lowg)), - col = col, border = NA) - points(x$r,x$theo,type=type,col=2) - points(x$r,x$aveg,type=type) -} From 258394a143221604eef8b4e418599745ff6f7d2c Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 17 Nov 2017 12:47:53 +0100 Subject: [PATCH 02/14] Trackstat.R --- R/Trackstat.R | 412 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 412 insertions(+) create mode 100644 R/Trackstat.R diff --git a/R/Trackstat.R b/R/Trackstat.R new file mode 100644 index 0000000..7f1010d --- /dev/null +++ b/R/Trackstat.R @@ -0,0 +1,412 @@ +# function as.Track accepts a dataframe (let say X) with three columns "xcoor", "ycoor" and time (or any class that can be converted to a data.frame) +# and converts it to an object of class Track. It can also accepts covariates for the corresponding +# locations, covariates must be a dataframe with some columns and length of each column is equal +# to length of each column in X. +as.Track <- function(X,covariate){ + stopifnot(nrow(X)>0) + colnames(X) <- c("xcoor","ycoor","time") + if(!is.data.frame(X)) X <- as.data.frame(X) + sp <- cbind(x=X$xcoor,y=X$ycoor) + sp <- SpatialPoints(sp) + t <- as.POSIXct(paste(X$date,X$time)) + if(missing(covariate)) covariate <- data.frame(d=rep(NA,length(X$xcoor))) + + return(Track(STIDF(sp,time = t,data =covariate))) +} + + +# function reTrack accepts X as an object of class Track. Output is a reconstructed Track (an object of class Track), based on "timestamp". +# It only returns the interpolated points. +reTrack <- function(X,at=c("track","dfrm"),timestamp=timestamp,tsq=NULL){ + + if (missing(tsq)) tsq <- tsqTracks(X,timestamp = timestamp) + if(missing(at)) at <- "track" + Xrange <- rngTrack(X) + X <- cbind(as.data.frame(X)[c(coordnames(X), "time")]) + xnew <- c() + ynew <- c() + + x <- X$x + y <- X$y + + time <- tsq[tsqXrange[1]] + ivs <- findInterval(time,X$time) + + for (i in 1:length(ivs)) { + if (!ivs[i] == 0 && !ivs[i] == nrow(X)) { + + iv <- ivs[i] + tdiff1 <- difftime(time[i],X$time[iv],units = "sec") # diff between timestamp and start of the interval it falls in + tdiff2 <- difftime(X$time[iv+1],X$time[iv],units = units(tdiff1)) # diff between timestamps (calculated here because it often varies) + ratio <- as.numeric(tdiff1)/as.numeric(tdiff2) + x1 <- X[iv,1] # segment coordinates + y1 <- X[iv,2] + x2 <- X[iv+1,1] + y2 <- X[iv+1,2] + xnew <- c(xnew, x1 + ratio * (x2 - x1)) #find point + ynew <- c(ynew, y1 + ratio * (y2 - y1)) + } + } + + newTrack <- data.frame(xnew, ynew, time) + newTrack <- newTrack[!duplicated(newTrack),] # remove duplicates + newTrack <- newTrack[order(newTrack$time),] # sort by timestamp + colnames(newTrack) <- c("xcoor","ycoor","time") + if (at=="dfrm") {attr(newTrack,"tsq") <-tsq;return(newTrack) } + return(as.Track(newTrack)) +} + +# rngTrack returns the timerange of an object of class Track +rngTrack <- function(X){ + Y <- cbind(as.data.frame(X)[c(coordnames(X), "time")]) + return(range(Y$time)) +} + +# tsqtracks returns a sequance of time based on a list of tracks (or a single object of class Track) and an argument timestamp +tsqTracks <- function(X,timestamp){ + + if (!is.list(X)) timerange <- rngTrack(X) + else timerange <- lapply(X, rngTrack) + + Trackrg <- range(timerange) + class(Trackrg) <- c('POSIXt','POSIXct') + # a seq from the range has been created every timestamp + timeseq <- seq(from=as.POSIXct(strftime(Trackrg[1])),to=as.POSIXct(strftime(Trackrg[2])),by = timestamp) + + return(timeseq) + +} + + + +# function avedistTrack accepts X as a list of tracks and reports the average distance between +# tracks over time, output is an object of class "distrack" +avedistTrack <- function(X,timestamp){ + + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + # calculate a sequance of time to interpolate tracks within this sequance + timeseq <- tsqTracks(X,timestamp = timestamp) + + # reconstruct tracks in sequance timeseq + Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") + + wincor <- lapply(X=1:length(Z),FUN = function(i){ + return(list(min(Z[[i]]$x),max(Z[[i]]$x),min(Z[[i]]$y),max(Z[[i]]$y))) + }) + wincor <- matrix(unlist(wincor),nrow = 4) + w <- owin(c(min(wincor[1,]),max(wincor[2,])),c(min(wincor[3,]),max(wincor[4,]))) + # create a list and convert tracks in each element of timeseq to an object of calss ppp + p <- list() + + for (j in 1:length(timeseq)) { + + l <- lapply(X=1:length(Z), function(i){ + Z[[i]][which(Z[[i]]$time==timeseq[j]),-3] + }) + + if(length(unlist(l))>0){ + x <- unlist(lapply(X=1:length(l),function(k){ + l[[k]]$x + })) + y <- unlist(lapply(X=1:length(l),function(h){ + l[[h]]$y + })) + p[[j]] <- as.ppp(data.frame(x,y),W=w) + } + + } + + p1 <- p[!sapply(p, is.null)] + avedist <- unlist( + lapply(X=1:length(p1), function(i){ + pd <- pairdist(p1[[i]]) + pd <- pd[pd>0] + return(mean(pd)) + }) + ) + avedist <- data.frame(timeseq[!sapply(p, is.null)],avedist) + class(avedist) <- c("distrack") + attr(avedist,"ppp") <- p + return(avedist) +} + +print.distrack <- function(x){ + print(as.vector(x$avedist)) +} + +plot.distrack <- function(x,...){ + plot(x$timeseq,x$avedist,,xlab="time",ylab="average distance",...) +} + + +rmdupTrack <- function(X){ + X <- cbind(as.data.frame(X)[c(coordnames(X), "time")]) + X <- X[!duplicated(X),] + return(as.Track(X)) +} + +as.Track.ppp <- function(X,timestamp){ + + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + # calculate a sequance of time to interpolate tracks within this sequance + timeseq <- tsqTracks(X,timestamp = timestamp) + + # reconstruct tracks in sequance timeseq + Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") + + wincor <- lapply(X=1:length(Z),FUN = function(i){ + return(list(min(Z[[i]]$x),max(Z[[i]]$x),min(Z[[i]]$y),max(Z[[i]]$y))) + }) + wincor <- matrix(unlist(wincor),nrow = 4) + w <- owin(c(min(wincor[1,]),max(wincor[2,])),c(min(wincor[3,]),max(wincor[4,]))) + + # create a list and convert tracks in each element of timeseq to an object of calss ppp + p <- list() + + for (j in 1:length(timeseq)) { + + l <- lapply(X=1:length(Z), function(i){ + w <- which(Z[[i]]$time==timeseq[j]) + if (length(w)>0) return(cbind(Z[[i]][w,-3],id=i)) + return(Z[[i]][w,-3]) + }) + + if(length(unlist(l))>0){ + x <- unlist(lapply(X=1:length(l),function(k){ + l[[k]]$x + })) + y <- unlist(lapply(X=1:length(l),function(h){ + l[[h]]$y + })) + m <- unlist(lapply(X=1:length(l),function(h){ + l[[h]]$id + })) + p[[j]] <- as.ppp(data.frame(x,y),W=w) + marks(p[[j]]) <- m + } + + } + return(p) +} + + +density.Track <- function(X,timestamp,...){ + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + + p <- as.Track.ppp(X,timestamp) + p <- p[!sapply(p, is.null)] + imlist <- lapply(p, density.ppp,...) + out <- Reduce("+",imlist)/length(imlist) + attr(out,"Tracksim") <- imlist + attr(out,"ppps") <- p + return(out) +} + +as.Track.arrow <- function(X,timestamp,epsilon=epsilon){ + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + + Z <- as.Track.ppp(X,timestamp) + Z <- Z[!sapply(Z, is.null)] + wind <- Z[[1]]$window + arrows <- list() + Y <- list() + for (i in 1:length(Z)) { + if(i==length(Z)) break() + j <- i+1 + m1 <- match(marks(Z[[i]]),marks(Z[[j]])) + m2 <- match(marks(Z[[j]]),marks(Z[[i]])) + m1 <- m1[!is.na(m1)] + m2 <- m2[!is.na(m2)] + x <- Z[[j]][m1] + y <- Z[[i]][m2] + l <- psp(y$x,y$y,x$x,x$y,window = wind) + arrows[[i]] <- l + center <- midpoints.psp(l) + mark <- lengths.psp(l) + marks(center) <- mark + if (missing(epsilon)) epsilon <- 0 + Y[[i]] <- center[mark>epsilon] + } + class(Y) <- c("list","Trrow") + attr(Y,"psp") <- arrows + return(Y) +} + +print.Trrow <- function(x, ...) { + attributes(x) <- NULL + print(x) +} + +Track.idw <- function(X,timestamp,epsilon=epsilon,...){ + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + + Y <- as.Track.arrow(X,timestamp,epsilon=epsilon) + Z <- lapply(Y, idw,...) + meanIDW <- Reduce("+",Z)/length(Z) + return(meanIDW) +} + +avemove <- function(X,timestamp,epsilon=epsilon){ + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + timeseq <- tsqTracks(X,timestamp = timestamp) + Y <- as.Track.arrow(X,timestamp,epsilon=epsilon) + Z <- attr(Y,"psp") + preout <- lapply(X=1:length(Z), function(i){ + mean(lengths.psp(Z[[i]])) + }) + out <- unlist(preout) + class(out) <- c("numeric", "arwlen") + attr(out,"tsq") <- timeseq + return(out) +} + +print.arwlen <- function(x){ + print(as.vector(x)) +} + +plot.arwlen <- function(x,...){ + tsq <- attr(x,"tsq") + tsq <- tsq[c(-1,-length(tsq))] + plot(tsq,x,xlab="time",ylab="average movement",...) +} + +chimaps <- function(X,timestamp,rank,...){ + if (!is.numeric(rank)) stop("rank must be numeric") + if (rank < 1 | rank >length(X)) stop("rank must be number between one and the number of Tracks") + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + timeseq <- tsqTracks(X,timestamp = timestamp) + d <- density.Track(X,timestamp,...) + imlist <- attr(d,"Tracksim") + sumim <- Reduce("+",imlist) + chi <- lapply(X=1:length(imlist),FUN = function(i){ + E1 <- sumim*sum(imlist[[i]]$v)/(sum(sumim$v)) + return((imlist[[i]]-E1)/sqrt(E1)) + }) + out <- chi[[rank]] + attr(out,"ims") <- chi + attr(out,"time") <- timeseq[rank] + attr(out,"timevec") <- timeseq + return(out) +} + +Kinhom.Track <- function(X,timestamp, + correction=c("border", "bord.modif", "isotropic", "translate"),q,...){ + + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + + cor <- match.arg(correction,correction) + + ZZ <- density.Track(X,timestamp) + + Z <- attr(ZZ,"Tracksim") + Y <- attr(ZZ,"ppps") + W <- Y[[1]]$window + ripley <- min(diff(W$xrange), diff(W$yrange))/4 + rr <- seq(0,ripley,length.out = 513) + + K <- lapply(X=1:length(Y), function(i){ + kk <- Kinhom(Y[[i]],lambda = Z[[i]],correction=cor,r=rr,...) + return(as.data.frame(kk)) + }) + Kmat <- matrix(nrow = length(K[[1]]$theo),ncol = length(K)) + for (i in 1:length(K)) { + Kmat[,i] <- K[[i]][,3] + } + # Kmat <- as.data.frame(K) + lowk <- numeric() + upk <- numeric() + avek <- numeric() + for (i in 1:nrow(Kmat)) { + avek[i] <- mean(Kmat[i,]) + lowk[i] <- quantile(Kmat[i,],q) + upk[i] <- quantile(Kmat[i,],1-q) + } + + out <- data.frame(lowk=lowk,upk=upk,avek=avek,r=K[[1]]$r,theo=K[[1]]$theo) + class(out) <- c("list","KTrack") + attr(out,"out") <- out + return(out) +} +print.KTrack <- function(x){ + print("variability area of K-function") +} + +plot.KTrack <- function(x,type="l",col= "grey70",...){ + ylim <- c(min(x$lowk),max(x$upk)) + plot(x$r,x$lowk,ylim=ylim,xlab="r",ylab="K(r)",type=type,...) + points(x$r,x$upk,type=type) + polygon(c(x$r, rev(x$r)), c(x$upk, rev(x$lowk)), + col = col, border = NA) + points(x$r,x$theo,type=type,col=2) + points(x$r,x$avek,type=type) + legend(0,max(x$upk),col = c(2,1),legend=c("poisson","average"),lty=c(1,1)) +} + +pcfinhom.Track <- function(X,timestamp, + correction = c("translate", "Ripley"),q,...){ + + stopifnot(length(X)>1 & is.list(X)) + + if (missing(timestamp)) stop("set timestamp") + + cor <- match.arg(correction,correction) + + ZZ <- density.Track(X,timestamp) + + Z <- attr(ZZ,"Tracksim") + Y <- attr(ZZ,"ppps") + + g <- lapply(X=1:length(Y), function(i){ + gg <- pcfinhom(Y[[i]],lambda = Z[[i]],correction=cor,...) + return(as.data.frame(gg)) + }) + gmat <- matrix(nrow = length(g[[1]]$theo),ncol = length(g)) + for (i in 1:length(g)) { + gmat[,i] <- g[[i]][,3] + } + # Kmat <- as.data.frame(K) + lowg <- numeric() + upg <- numeric() + aveg <- numeric() + for (i in 1:nrow(gmat)) { + aveg[i] <- mean(gmat[i,]) + lowg[i] <- quantile(gmat[i,],q) + upg[i] <- quantile(gmat[i,],1-q) + } + + out <- data.frame(lowg=lowg,upg=upg,aveg=aveg,r=g[[1]]$r,theo=g[[1]]$theo) + class(out) <- c("list","gTrack") + attr(out,"out") <- out + return(out) +} + + +print.gTrack <- function(x){ + print("variability area of pair correlatio function") +} + +plot.gTrack <- function(x,type="l",col= "grey70",...){ + ylim <- c(min(x$lowg),max(x$upg)) + plot(x$r,x$lowg,ylim=ylim,xlab="r",ylab="g",type=type,...) + points(x$r,x$upg,type=type) + polygon(c(x$r, rev(x$r)), c(x$upg, rev(x$lowg)), + col = col, border = NA) + points(x$r,x$theo,type=type,col=2) + points(x$r,x$aveg,type=type) +} From db284b24b72ae2cba95676867e7000d2711d4c62 Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 17 Nov 2017 17:00:11 +0100 Subject: [PATCH 03/14] improved as.Track.ppp --- R/Trackstat.R | 89 +++++++++++---------------------------------------- 1 file changed, 19 insertions(+), 70 deletions(-) diff --git a/R/Trackstat.R b/R/Trackstat.R index 7f1010d..29a0028 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -89,49 +89,19 @@ avedistTrack <- function(X,timestamp){ # calculate a sequance of time to interpolate tracks within this sequance timeseq <- tsqTracks(X,timestamp = timestamp) - # reconstruct tracks in sequance timeseq - Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") + Y <- as.Track.ppp(X,timestamp) - wincor <- lapply(X=1:length(Z),FUN = function(i){ - return(list(min(Z[[i]]$x),max(Z[[i]]$x),min(Z[[i]]$y),max(Z[[i]]$y))) + avedist <- lapply(X=1:length(Y), function(i){ + pd <- pairdist(Y[[i]]) + mean(pd[pd>0]) }) - wincor <- matrix(unlist(wincor),nrow = 4) - w <- owin(c(min(wincor[1,]),max(wincor[2,])),c(min(wincor[3,]),max(wincor[4,]))) - # create a list and convert tracks in each element of timeseq to an object of calss ppp - p <- list() - - for (j in 1:length(timeseq)) { - - l <- lapply(X=1:length(Z), function(i){ - Z[[i]][which(Z[[i]]$time==timeseq[j]),-3] - }) - - if(length(unlist(l))>0){ - x <- unlist(lapply(X=1:length(l),function(k){ - l[[k]]$x - })) - y <- unlist(lapply(X=1:length(l),function(h){ - l[[h]]$y - })) - p[[j]] <- as.ppp(data.frame(x,y),W=w) - } - - } - p1 <- p[!sapply(p, is.null)] - avedist <- unlist( - lapply(X=1:length(p1), function(i){ - pd <- pairdist(p1[[i]]) - pd <- pd[pd>0] - return(mean(pd)) - }) - ) - avedist <- data.frame(timeseq[!sapply(p, is.null)],avedist) + avedist <- data.frame(timeseq[-1],unlist(avedist)) + colnames(avedist) <- c("timeseq","avedist") class(avedist) <- c("distrack") - attr(avedist,"ppp") <- p + attr(avedist,"ppp") <- Y return(avedist) } - print.distrack <- function(x){ print(as.vector(x$avedist)) } @@ -158,42 +128,21 @@ as.Track.ppp <- function(X,timestamp){ # reconstruct tracks in sequance timeseq Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") - wincor <- lapply(X=1:length(Z),FUN = function(i){ - return(list(min(Z[[i]]$x),max(Z[[i]]$x),min(Z[[i]]$y),max(Z[[i]]$y))) + Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") + id <- rep(1:length(Z),sapply(Z, nrow)) + Z <- do.call("rbind",Z) + Z <- cbind(Z,id) + allZ <- split(Z,Z[,3]) + w <- owin(c(min(Z$xcoor)-0.5,max(Z$xcoor)+0.5),c(min(Z$ycoor)-0.5,max(Z$ycoor)+0.5)) + + Tppp <- lapply(X=1:length(allZ), function(i){ + p <- as.ppp(allZ[[i]][,-c(3,4)],W=w) + marks(p) <- allZ[[i]][,4] + return(p) }) - wincor <- matrix(unlist(wincor),nrow = 4) - w <- owin(c(min(wincor[1,]),max(wincor[2,])),c(min(wincor[3,]),max(wincor[4,]))) - - # create a list and convert tracks in each element of timeseq to an object of calss ppp - p <- list() - - for (j in 1:length(timeseq)) { - - l <- lapply(X=1:length(Z), function(i){ - w <- which(Z[[i]]$time==timeseq[j]) - if (length(w)>0) return(cbind(Z[[i]][w,-3],id=i)) - return(Z[[i]][w,-3]) - }) - - if(length(unlist(l))>0){ - x <- unlist(lapply(X=1:length(l),function(k){ - l[[k]]$x - })) - y <- unlist(lapply(X=1:length(l),function(h){ - l[[h]]$y - })) - m <- unlist(lapply(X=1:length(l),function(h){ - l[[h]]$id - })) - p[[j]] <- as.ppp(data.frame(x,y),W=w) - marks(p[[j]]) <- m - } - - } - return(p) + return(Tppp) } - density.Track <- function(X,timestamp,...){ stopifnot(length(X)>1 & is.list(X)) From 9205c092379fd97b110ae8e8cdec5df9aee03f06 Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 17 Nov 2017 17:45:16 +0100 Subject: [PATCH 04/14] minor edit --- R/Trackstat.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/Trackstat.R b/R/Trackstat.R index 29a0028..b326646 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -126,8 +126,6 @@ as.Track.ppp <- function(X,timestamp){ timeseq <- tsqTracks(X,timestamp = timestamp) # reconstruct tracks in sequance timeseq - Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") - Z <- lapply(X,reTrack,tsq = timeseq,at="dfrm") id <- rep(1:length(Z),sapply(Z, nrow)) Z <- do.call("rbind",Z) From 2bf6dc2538765ee0a9536088177141bf7c4ea571 Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 17 Nov 2017 18:13:01 +0100 Subject: [PATCH 05/14] edit description --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index cde0b02..258e726 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -17,4 +17,4 @@ URL: http://github.com/edzer/trajectories BugReports: http://github.com/edzer/trajectories/issues VignetteBuilder: knitr Collate: Class-Tracks.R Tracks-methods.R generalize.R stcube.R stplot.R - difftrack.R compare-methods.R rtracks.R + difftrack.R compare-methods.R rtracks.R Trackstat.R From 1abdd7ab8d7a99b187c48dcabde5b30d1a4db8e6 Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 17 Nov 2017 18:13:24 +0100 Subject: [PATCH 06/14] to check --- R/Trackstat.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/Trackstat.R b/R/Trackstat.R index b326646..755b2b7 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -117,6 +117,7 @@ rmdupTrack <- function(X){ return(as.Track(X)) } + as.Track.ppp <- function(X,timestamp){ stopifnot(length(X)>1 & is.list(X)) From 7a21d78feff9c290bb03bce24040aaa41b4a5c04 Mon Sep 17 00:00:00 2001 From: Moradii Date: Thu, 30 Nov 2017 15:02:28 +0100 Subject: [PATCH 07/14] bug fixed --- R/Trackstat.R | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/R/Trackstat.R b/R/Trackstat.R index 755b2b7..b43ca9e 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -9,7 +9,7 @@ as.Track <- function(X,covariate){ sp <- cbind(x=X$xcoor,y=X$ycoor) sp <- SpatialPoints(sp) t <- as.POSIXct(paste(X$date,X$time)) - if(missing(covariate)) covariate <- data.frame(d=rep(NA,length(X$xcoor))) + if(missing(covariate)) covariate <- data.frame(d=rep(0,length(X$xcoor))) return(Track(STIDF(sp,time = t,data =covariate))) } @@ -132,7 +132,7 @@ as.Track.ppp <- function(X,timestamp){ Z <- do.call("rbind",Z) Z <- cbind(Z,id) allZ <- split(Z,Z[,3]) - w <- owin(c(min(Z$xcoor)-0.5,max(Z$xcoor)+0.5),c(min(Z$ycoor)-0.5,max(Z$ycoor)+0.5)) + w <- owin(c(min(Z$xcoor)-0.001,max(Z$xcoor)+0.001),c(min(Z$ycoor)-0.001,max(Z$ycoor)+0.001)) Tppp <- lapply(X=1:length(allZ), function(i){ p <- as.ppp(allZ[[i]][,-c(3,4)],W=w) @@ -252,15 +252,17 @@ chimaps <- function(X,timestamp,rank,...){ } Kinhom.Track <- function(X,timestamp, - correction=c("border", "bord.modif", "isotropic", "translate"),q,...){ + correction=c("border", "bord.modif", "isotropic", "translate"),q, + sigma=c("bw.diggle","bw.ppl"," bw.scott"),...){ stopifnot(length(X)>1 & is.list(X)) if (missing(timestamp)) stop("set timestamp") cor <- match.arg(correction,correction) - - ZZ <- density.Track(X,timestamp) + bw <- match.arg(sigma,sigma) + bw <- match.fun(bw) + ZZ <- density.Track(X,timestamp,bw) Z <- attr(ZZ,"Tracksim") Y <- attr(ZZ,"ppps") From 5e70c8cacf4daadb1b04b2f3099aba6691786f32 Mon Sep 17 00:00:00 2001 From: Moradii Date: Wed, 6 Dec 2017 15:29:45 +0100 Subject: [PATCH 08/14] newrTrack,newrTracks,newrTracksCollection these three functions accept box for simulating tracks --- R/Trackstat.R | 52 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/R/Trackstat.R b/R/Trackstat.R index b43ca9e..83cf391 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -360,3 +360,55 @@ plot.gTrack <- function(x,type="l",col= "grey70",...){ points(x$r,x$theo,type=type,col=2) points(x$r,x$aveg,type=type) } + + +newrTrack <- function (n = 100, origin = c(0, 0), start = as.POSIXct("1970-01-01"), + ar = 0.8, step = 60, sd0 = 1,bbox=bbox, transform=FALSE, ...){ + if (missing(bbox) & transform) { + xo <- runif(1) + yo <- runif(1) + origin <- c(xo,yo) + } + if (!missing(bbox) & transform) { + xo <- runif(1,bbox[1,1],bbox[1,2]) + yo <- runif(1,bbox[2,1],bbox[2,2]) + origin <- c(xo,yo) + } + if (length(ar) == 1 && ar == 0) + xy = cbind(cumsum(rnorm(n, sd = sd0)) + origin[1], cumsum(rnorm(n, + sd = sd0)) + origin[2]) + else {xy = cbind(origin[1] + cumsum(as.vector(arima.sim(list(ar = ar), + n, sd = sd0, ...))), + origin[2] + cumsum(as.vector(arima.sim(list(ar = ar), + + n, sd = sd0, ...))))} + if(transform) { + if(missing(bbox)) bbox <- matrix(c(0,1,0,1),nrow = 2,byrow = T); colnames(bbox) <- c("min","max");rownames(bbox) <- c("x","y") + + xr <- max(xy[,1])-min(xy[,1]) + yr <- max(xy[,2])-min(xy[,2]) + + xt <- (xy[,1]-min(xy[,1]))/xr + yt <- (xy[,2]-min(xy[,2]))/yr + + xy <- cbind(xt,yt) + xy <- cbind(x=xy[,1]*bbox[1,2],y=xy[,2]*bbox[2,2]) + } + + T = start + 0:(n - 1) * step + sti = STI(SpatialPoints(xy), T) + out <- Track(sti) + if (transform) out@sp@bbox <- bbox + return(out) +} + + +newrTracks <- function (m = 20, start = as.POSIXct("1970-01-01"), delta = 7200, + sd1 = 0, origin = c(0, 0), ...) + Tracks(lapply(0:(m - 1) * delta, function(x) newrTrack(start = start + + x, origin = origin + rnorm(2, sd = sd1), ...))) + +newrTracksCollection <- function (p = 10, sd2 = 0, ...) + TracksCollection(lapply(1:p, function(x) newrTracks(origin = rnorm(2, + sd = sd2), ...))) + From 105aadd7ef1fd2446013185cd6646da5216489b0 Mon Sep 17 00:00:00 2001 From: Moradii Date: Thu, 7 Dec 2017 14:08:07 +0100 Subject: [PATCH 09/14] arg "nrandom" is added this allows the number of points per track to be a random number from Poisson dist. --- R/Trackstat.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/Trackstat.R b/R/Trackstat.R index 83cf391..13382ab 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -363,7 +363,9 @@ plot.gTrack <- function(x,type="l",col= "grey70",...){ newrTrack <- function (n = 100, origin = c(0, 0), start = as.POSIXct("1970-01-01"), - ar = 0.8, step = 60, sd0 = 1,bbox=bbox, transform=FALSE, ...){ + ar = 0.8, step = 60, sd0 = 1,bbox=bbox, transform=FALSE,nrandom=FALSE, ...){ + + if(nrandom) repeat{n <- rpois(1,n);if(!a==0) break()} if (missing(bbox) & transform) { xo <- runif(1) yo <- runif(1) From 3de48ad76edb67193d36eb48b8e6d014d6f952a5 Mon Sep 17 00:00:00 2001 From: Moradii Date: Thu, 7 Dec 2017 14:41:25 +0100 Subject: [PATCH 10/14] bug fixed --- R/Trackstat.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Trackstat.R b/R/Trackstat.R index 13382ab..0fd6652 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -365,7 +365,7 @@ plot.gTrack <- function(x,type="l",col= "grey70",...){ newrTrack <- function (n = 100, origin = c(0, 0), start = as.POSIXct("1970-01-01"), ar = 0.8, step = 60, sd0 = 1,bbox=bbox, transform=FALSE,nrandom=FALSE, ...){ - if(nrandom) repeat{n <- rpois(1,n);if(!a==0) break()} + if(nrandom) repeat{n <- rpois(1,n);if(!n==0) break()} if (missing(bbox) & transform) { xo <- runif(1) yo <- runif(1) From 903cf5ac6af3b31a36c42b55618f137efd2eae63 Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 8 Dec 2017 17:31:23 +0100 Subject: [PATCH 11/14] keeping rtracks with the pervious name --- R/Trackstat.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/Trackstat.R b/R/Trackstat.R index 0fd6652..402b3d7 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -4,12 +4,12 @@ # to length of each column in X. as.Track <- function(X,covariate){ stopifnot(nrow(X)>0) - colnames(X) <- c("xcoor","ycoor","time") + # colnames(X) <- c("xcoor","ycoor","time") if(!is.data.frame(X)) X <- as.data.frame(X) sp <- cbind(x=X$xcoor,y=X$ycoor) sp <- SpatialPoints(sp) t <- as.POSIXct(paste(X$date,X$time)) - if(missing(covariate)) covariate <- data.frame(d=rep(0,length(X$xcoor))) + if(missing(covariate)) covariate <- data.frame(d=rep(1,length(X$xcoor))) return(Track(STIDF(sp,time = t,data =covariate))) } @@ -362,7 +362,7 @@ plot.gTrack <- function(x,type="l",col= "grey70",...){ } -newrTrack <- function (n = 100, origin = c(0, 0), start = as.POSIXct("1970-01-01"), +rTrack <- function (n = 100, origin = c(0, 0), start = as.POSIXct("1970-01-01"), ar = 0.8, step = 60, sd0 = 1,bbox=bbox, transform=FALSE,nrandom=FALSE, ...){ if(nrandom) repeat{n <- rpois(1,n);if(!n==0) break()} @@ -405,12 +405,12 @@ newrTrack <- function (n = 100, origin = c(0, 0), start = as.POSIXct("1970-01-01 } -newrTracks <- function (m = 20, start = as.POSIXct("1970-01-01"), delta = 7200, +rTracks <- function (m = 20, start = as.POSIXct("1970-01-01"), delta = 7200, sd1 = 0, origin = c(0, 0), ...) Tracks(lapply(0:(m - 1) * delta, function(x) newrTrack(start = start + x, origin = origin + rnorm(2, sd = sd1), ...))) -newrTracksCollection <- function (p = 10, sd2 = 0, ...) +rTracksCollection <- function (p = 10, sd2 = 0, ...) TracksCollection(lapply(1:p, function(x) newrTracks(origin = rnorm(2, sd = sd2), ...))) From f324a479ce0c03d1471af1fc856acb4b4621ba4c Mon Sep 17 00:00:00 2001 From: Moradii Date: Fri, 8 Dec 2017 19:57:01 +0100 Subject: [PATCH 12/14] bug fixed --- R/Trackstat.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/Trackstat.R b/R/Trackstat.R index 402b3d7..5222e77 100644 --- a/R/Trackstat.R +++ b/R/Trackstat.R @@ -407,10 +407,10 @@ rTrack <- function (n = 100, origin = c(0, 0), start = as.POSIXct("1970-01-01"), rTracks <- function (m = 20, start = as.POSIXct("1970-01-01"), delta = 7200, sd1 = 0, origin = c(0, 0), ...) - Tracks(lapply(0:(m - 1) * delta, function(x) newrTrack(start = start + + Tracks(lapply(0:(m - 1) * delta, function(x) rTrack(start = start + x, origin = origin + rnorm(2, sd = sd1), ...))) rTracksCollection <- function (p = 10, sd2 = 0, ...) - TracksCollection(lapply(1:p, function(x) newrTracks(origin = rnorm(2, + TracksCollection(lapply(1:p, function(x) rTracks(origin = rnorm(2, sd = sd2), ...))) From 0edfe94658486433ccd31748e2ed7fdba58d2649 Mon Sep 17 00:00:00 2001 From: Edzer Pebesma Date: Mon, 11 Dec 2017 14:30:42 +0100 Subject: [PATCH 13/14] fix docs and namespace issues --- DESCRIPTION | 2 +- NAMESPACE | 8 ++++++-- man/Track-class.Rd | 3 +++ man/rtrack.Rd | 5 ++++- 4 files changed, 14 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 258e726..54b04a9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,7 +8,7 @@ Authors@R: person("Benedikt", "Graeler", role = "ctb"), person("Nikolai", "Gorte", role = "ctb")) Depends: R (>= 3.0.0) -Imports: stats, utils, graphics, methods, lattice, sp (>= 1.1-0), spacetime (>= 1.0-0) +Imports: stats, utils, graphics, methods, lattice, sp (>= 1.1-0), spacetime (>= 1.0-0), spatstat Suggests: rgdal, rgeos, rgl, OpenStreetMap, RCurl, rjson, adehabitatLT, xts, knitr LazyData: no Description: Classes and methods for trajectory data, with support for nesting individual Track objects in track sets (Tracks) and track sets for different entities in collections of Tracks (TracksCollection). Methods include selection, generalization, aggregation, intersection, simulation, and plotting. diff --git a/NAMESPACE b/NAMESPACE index fdf1ea7..8de8eca 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,11 +2,15 @@ import(methods) import(sp) import(spacetime) import(lattice) -importFrom(stats, approx, as.formula, na.omit, quantile, rnorm, arima.sim) +importFrom(stats, approx, as.formula, na.omit, quantile, rnorm, arima.sim, + rpois, runif) # importFrom(stats, aggregate, na.omit, time, start, end) importFrom(utils, head, tail, stack, unstack) -importFrom(graphics, arrows, box, lines, points, segments) +importFrom(graphics, arrows, box, lines, points, segments, legend, polygon) importFrom(grDevices, rainbow) +importFrom(spatstat, marks, psp, midpoints.psp, lengths.psp, + 'marks<-', owin, as.ppp, pairdist, density.psp, density.ppp, + Kinhom, pcfinhom, idw) exportClasses( Track, diff --git a/man/Track-class.Rd b/man/Track-class.Rd index c04e6a8..216ef11 100644 --- a/man/Track-class.Rd +++ b/man/Track-class.Rd @@ -13,6 +13,9 @@ \alias{[,Track-method} \alias{[,Tracks-method} \alias{[,TracksCollection-method} +\alias{[,Track,ANY,ANY,ANY-method} +\alias{[,Tracks,ANY,ANY,ANY-method} +\alias{[,TracksCollection,ANY,ANY,ANY-method} \alias{[[,Track,ANY,missing-method} \alias{[[,Tracks,ANY,missing-method} \alias{[[,TracksCollection,ANY,missing-method} diff --git a/man/rtrack.Rd b/man/rtrack.Rd index db46529..bfb25ad 100644 --- a/man/rtrack.Rd +++ b/man/rtrack.Rd @@ -8,7 +8,7 @@ \description{Generate random \code{Track}, \code{Tracks} or \code{TracksCollection} objects} \usage{ rTrack(n = 100, origin = c(0,0), start = as.POSIXct("1970-01-01"), ar = .8, - step = 60, sd0 = 1, ...) + step = 60, sd0 = 1, bbox = bbox, transform = FALSE, nrandom = FALSE, ...) rTracks(m = 20, start = as.POSIXct("1970-01-01"), delta = 7200, sd1 = 0, origin = c(0,0), ...) rTracksCollection(p = 10, sd2 = 0, ...) @@ -22,6 +22,9 @@ rTracksCollection(p = 10, sd2 = 0, ...) \item{sd0}{standard deviation of the random steps in a Track} \item{sd1}{standard deviation of the consecutive Track origin values (using rnorm)} \item{sd2}{standard deviation of the consecutive Tracks origin values (using rnorm)} +\item{bbox}{bbox object FIXME:fill in} +\item{transform}{logical; FIXME:fill in } +\item{nrandom}{logical; if \code{TRUE}, draw \code{n} from \code{rpois(n)}} \item{...}{rTrack: arguments passed on to \link[stats]{arima.sim}, rTracks: arguments passed on to rTrack; rTracksCollection: arguments passed on to rTracks} \item{m}{ number of Track objects to simulate} From 7c862e2d3469e2ceaaa75321fb80f05f90f18aa1 Mon Sep 17 00:00:00 2001 From: Edzer Pebesma Date: Thu, 29 Mar 2018 21:41:24 +0200 Subject: [PATCH 14/14] add vignette --- vignettes/trsim.html | 389 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 389 insertions(+) create mode 100644 vignettes/trsim.html diff --git a/vignettes/trsim.html b/vignettes/trsim.html new file mode 100644 index 0000000..8cde039 --- /dev/null +++ b/vignettes/trsim.html @@ -0,0 +1,389 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + + +
+

Interpolating, smoothing, and simulating trajectories

+
+

create some data

+
t0 = as.POSIXct(as.Date("2013-09-30",tz="CET"))
+# person A, track 1:
+x = c(7,6,5,5,4,3,3)
+y = c(7,7,6,5,5,6,7)
+n = length(x)
+set.seed(131)
+t = t0 + cumsum(runif(n) * 60)
+library(sp)
+require(rgdal)
+
## Loading required package: rgdal
+
## rgdal: version: 1.2-15, (SVN revision 691)
+##  Geospatial Data Abstraction Library extensions to R successfully loaded
+##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
+##  Path to GDAL shared files: /usr/share/gdal/2.1
+##  GDAL binary built with GEOS: TRUE 
+##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
+##  Path to PROJ.4 shared files: (autodetected)
+##  Linking to sp version: 1.2-5
+
crs = CRS("+proj=longlat +ellps=WGS84") # longlat
+
+library(spacetime)
+stidf = STIDF(SpatialPoints(cbind(x,y),crs), t, data.frame(co2 = rnorm(n)))
+
+library(trajectories)
+A1 = Track(stidf)
+# person A, track 2:
+x = c(7,6,6,7,7)
+y = c(6,5,4,4,3)
+n = length(x)
+t = max(t) + cumsum(runif(n) * 60)
+stidf = STIDF(SpatialPoints(cbind(x,y),crs), t, data.frame(co2 = rnorm(n)))
+
+A2 = Track(stidf)
+# Tracks for person A:
+A = Tracks(list(A1=A1,A2=A2))
+# person B, track 1:
+x = c(2,2,1,1,2,3)
+y = c(5,4,3,2,2,3)
+n = length(x)
+t = max(t) + cumsum(runif(n) * 60)
+stidf = STIDF(SpatialPoints(cbind(x,y),crs), t, data.frame(co2 = rnorm(n)))
+B1 = Track(stidf)
+# person B, track 2:
+x = c(3,3,4,3,3,4)
+y = c(5,4,3,2,1,1)
+n = length(x)
+t = max(t) + cumsum(runif(n) * 60)
+stidf = STIDF(SpatialPoints(cbind(x,y),crs), t, data.frame(co2 = rnorm(n)))
+B2 = Track(stidf)
+# Tracks for person A:
+B = Tracks(list(B1=B1,B2=B2))
+Tr = TracksCollection(list(A=A,B=B))
+stplot(Tr, scales = list(draw=TRUE))
+

+
stplot(Tr, attr = "direction", arrows=TRUE, lwd = 3, by = "direction")
+

+
stplot(Tr, attr = "direction", arrows=TRUE, lwd = 3, by = "IDs")
+

+
plot(Tr, col=2, axes=TRUE)
+

+
dim(Tr)
+
##        IDs     tracks geometries 
+##          2          4         24
+
dim(Tr[2])
+
##     tracks geometries 
+##          2         12
+
dim(Tr[2][1])
+
## geometries 
+##          6
+
u = stack(Tr) # four IDs
+dim(u)
+
##        IDs     tracks geometries 
+##          4          4         24
+
dim(unstack(u, c(1,1,2,2))) # regroups to original
+
##        IDs     tracks geometries 
+##          2          4         24
+
dim(unstack(u, c(1,1,2,3))) # regroups to three IDs
+
##        IDs     tracks geometries 
+##          3          4         24
+
dim(unstack(u, c(1,2,2,1))) # regroups differently
+
##        IDs     tracks geometries 
+##          2          4         24
+
as(Tr, "data.frame")[1:10,] # tracks separated by NA rows
+
##         x  y sp.ID                time             endTime timeIndex
+## A.A1.1  7  7     1 2013-09-30 02:00:12 2013-09-30 02:00:12         1
+## A.A1.2  6  7     2 2013-09-30 02:00:19 2013-09-30 02:00:19         2
+## A.A1.3  5  6     3 2013-09-30 02:00:37 2013-09-30 02:00:37         3
+## A.A1.4  5  5     4 2013-09-30 02:01:00 2013-09-30 02:01:00         4
+## A.A1.5  4  5     5 2013-09-30 02:01:50 2013-09-30 02:01:50         5
+## A.A1.6  3  6     6 2013-09-30 02:02:22 2013-09-30 02:02:22         6
+## A.A1.7  3  7     7 2013-09-30 02:02:53 2013-09-30 02:02:53         7
+## A.A1.8 NA NA    NA                <NA>                <NA>        NA
+## A.A2.1  7  6     1 2013-09-30 02:03:31 2013-09-30 02:03:31         1
+## A.A2.2  6  5     2 2013-09-30 02:04:09 2013-09-30 02:04:09         2
+##                co2 Track IDs
+## A.A1.1 -0.71322105    A1   A
+## A.A1.2  1.37185185    A1   A
+## A.A1.3 -0.39982855    A1   A
+## A.A1.4 -0.47880016    A1   A
+## A.A1.5 -0.54870456    A1   A
+## A.A1.6  0.48757652    A1   A
+## A.A1.7 -0.06981164    A1   A
+## A.A1.8          NA    A1   A
+## A.A2.1  1.40377616    A2   A
+## A.A2.2  0.79969808    A2   A
+
as(Tr, "segments")[1:10,]   # track segments as records
+
##        x0 y0 x1 y1                time        co2 distance duration
+## A.A1.1  7  7  6  7 2013-09-30 02:00:12 -0.7132211 110495.2  7.49653
+## A.A1.2  6  7  5  6 2013-09-30 02:00:19  1.3718518 156407.3 17.59639
+## A.A1.3  5  6  5  5 2013-09-30 02:00:37 -0.3998285 110583.3 22.54678
+## A.A1.4  5  5  4  5 2013-09-30 02:01:00 -0.4788002 110898.7 50.78081
+## A.A1.5  4  5  3  6 2013-09-30 02:01:50 -0.5487046 156547.2 31.75229
+## A.A1.6  3  6  3  7 2013-09-30 02:02:22  0.4875765 110587.4 31.11752
+## A.A2.1  7  6  6  5 2013-09-30 02:03:31  1.4037762 156547.2 37.52483
+## A.A2.2  6  5  6  4 2013-09-30 02:04:09  0.7996981 110579.9 24.03699
+## A.A2.3  6  4  7  4 2013-09-30 02:04:33 -0.7110605 111050.1 34.96193
+## A.A2.4  7  4  7  3 2013-09-30 02:05:08 -1.5005960 110577.2 19.64212
+##            speed direction Track IDs
+## A.A1.1 14739.515 270.06094    A1   A
+## A.A1.2  8888.604 224.87295    A1   A
+## A.A1.3  4904.618 180.00000    A1   A
+## A.A1.4  2183.870 270.04358    A1   A
+## A.A1.5  4930.266 315.17903    A1   A
+## A.A1.6  3553.863   0.00000    A1   A
+## A.A2.1  4171.830 224.91682    A2   A
+## A.A2.2  4600.407 180.00000    A2   A
+## A.A2.3  3176.316  89.96512    A2   A
+## A.A2.4  5629.598 180.00000    A2   A
+
Tr[["distance"]] = Tr[["distance"]] * 1000
+Tr$distance = Tr$distance / 1000
+Tr$distance
+
##    A.A11    A.A12    A.A13    A.A14    A.A15    A.A16    A.A21    A.A22 
+## 110495.2 156407.3 110583.3 110898.7 156547.2 110587.4 156547.2 110579.9 
+##    A.A23    A.A24    B.B11    B.B12    B.B13    B.B14    B.B15    B.B21 
+## 111050.1 110577.2 110579.9 156757.4 110575.2 111252.1 156827.6 110579.9 
+##    B.B22    B.B23    B.B24    B.B25 
+## 156757.4 156827.6 110573.8 111302.6
+
# work with custum TrackStats function:
+MyStats = function(track) {
+    df = apply(coordinates(track@sp), 2, diff) # requires sp
+    data.frame(distance = apply(df, 1, function(x) sqrt(sum(x^2))))
+}
+crs = CRS(as.character(NA))
+stidf = STIDF(SpatialPoints(cbind(x,y),crs), t, data.frame(co2 = rnorm(n)))
+B2 = Track(stidf) # no longer longlat;
+B3 = Track(stidf, fn = MyStats)
+all.equal(B3$distance, B2$distance)
+
## [1] TRUE
+
+
+

Interpolating trajectories

+
opar = par()
+par(mfrow = c(1, 2))
+plot(B2, ylim = c(.5, 6))
+plot(B2, pch = 16, add = TRUE)
+title("irregular time steps")
+i = index(B2)
+B3 = approxTrack(B2, seq(min(i), max(i), length.out = 50))
+plot(B3, col = 'red', type = 'p', add = TRUE)
+B4 = approxTrack(B2, seq(min(i), max(i), length.out = 50), FUN = spline)
+plot(B4, col = 'blue', type = 'b', add = TRUE)
+# regular time steps:
+t = max(t) + (1:n) * 60 # regular
+B2 = Track(STIDF(SpatialPoints(cbind(x,y),crs), t, data.frame(co2 = rnorm(n))))
+plot(B2, ylim = c(.5, 6))
+plot(B2, pch = 16, add = TRUE)
+title("constant time steps")
+i = index(B2)
+B3 = approxTrack(B2)
+plot(B3, type = 'p', col = 'red', add = TRUE)
+B4 = approxTrack(B2, FUN = spline)
+plot(B4, type = 'p', col = 'blue', add = TRUE)
+

+
par(opar)
+
## Warning in par(opar): graphical parameter "cin" cannot be set
+
## Warning in par(opar): graphical parameter "cra" cannot be set
+
## Warning in par(opar): graphical parameter "csi" cannot be set
+
## Warning in par(opar): graphical parameter "cxy" cannot be set
+
## Warning in par(opar): graphical parameter "din" cannot be set
+
## Warning in par(opar): graphical parameter "page" cannot be set
+
+
+

Smoothing trajectories

+
smth = function(x,y,xout,...) predict(smooth.spline(as.numeric(x), y), as.numeric(xout))
+data(storms)
+plot(storms, type = 'p')
+
## Warning in axis(side, at = at, labels = labels, ...): graphical parameter
+## "type" is obsolete
+
+## Warning in axis(side, at = at, labels = labels, ...): graphical parameter
+## "type" is obsolete
+
## Warning in title(...): graphical parameter "type" is obsolete
+
storms.smooth = approxTracksCollection(storms, FUN = smth, n = 200)
+
## Warning in validityMethod(object): tracks with overlapping time intervals:
+## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18
+
## Warning in validityMethod(object): tracks with overlapping time intervals:
+## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
+
## Warning in validityMethod(object): tracks with overlapping time intervals:
+## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
+
## Warning in validityMethod(object): tracks with overlapping time intervals:
+## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
+
plot(storms.smooth, add = TRUE, col = 'red')
+

+
+
+

Simulating random trajectories

+
x = rTrack()
+dim(x)
+
## geometries 
+##        100
+
plot(x)
+

+
x = rTracks(sd1 = 120)
+dim(x)
+
##     tracks geometries 
+##         20       2000
+
plot(as(x, "SpatialLines"), col = 1:dim(x)[1], axes=TRUE)
+

+
x = rTracksCollection() # star
+dim(x)
+
##        IDs     tracks geometries 
+##         10        200      20000
+
plot(x)
+

+
x = rTracksCollection(sd2 = 200)
+plot(x, col=1:dim(x)[1])
+

+
+
+ + + + +
+ + + + + + + +