https://github.com/cran/spatstat
Tip revision: ad0cdcc7b4346a4492f3f16c48c26a8a6b710c8a authored by Adrian Baddeley on 09 May 2014, 11:24:31 UTC
version 1.37-0
version 1.37-0
Tip revision: ad0cdcc
crossdistlpp.R
#
# crossdistlpp.R
#
# $Revision: 1.3 $ $Date: 2013/08/23 07:35:36 $
#
# crossdist.lpp
# Calculates the shortest-path distance from each point of X
# to each point of Y, where X and Y are point patterns
# on the same linear network.
#
crossdist.lpp <- function(X, Y, ..., method="C") {
stopifnot(inherits(X, "lpp"))
stopifnot(method %in% c("C", "interpreted"))
check <- resolve.defaults(list(...), list(check=TRUE))$check
#
nX <- npoints(X)
nY <- npoints(Y)
#
L <- as.linnet(X)
if(check) {
LY <- as.linnet(Y)
if(!identical(L, LY))
stop("X and Y are on different linear networks")
}
P <- as.ppp(X)
Q <- as.ppp(Y)
#
Lseg <- L$lines
Lvert <- L$vertices
from <- L$from
to <- L$to
dpath <- L$dpath
# nearest segment for each point
Xpro <- coords(X, local=TRUE, spatial=FALSE, temporal=FALSE)$seg
Ypro <- coords(Y, local=TRUE, spatial=FALSE, temporal=FALSE)$seg
crossdistmat <- matrix(0,nX,nY)
if(method == "interpreted") {
# loop through all pairs of data points
for (i in 1:nX) {
Xproi <- Xpro[i]
Xi <- P[i]
nbi1 <- from[Xproi]
nbi2 <- to[Xproi]
vi1 <- Lvert[nbi1]
vi2 <- Lvert[nbi2]
dXi1 <- crossdist(Xi, vi1)
dXi2 <- crossdist(Xi, vi2)
for (j in 1:nY) {
Yj <- Q[j]
Yproj <- Ypro[j]
if(Xproi == Yproj) {
# points i and j lie on the same segment
# use Euclidean distance
d <- crossdist(Xi, Yj)
} else {
# shortest path from i to j passes through ends of segments
nbj1 <- from[Yproj]
nbj2 <- to[Yproj]
vj1 <- Lvert[nbj1]
vj2 <- Lvert[nbj2]
# Calculate shortest of 4 possible paths from i to j
d1Yj <- crossdist(vj1,Yj)
d2Yj <- crossdist(vj2,Yj)
d11 <- dXi1 + dpath[nbi1,nbj1] + d1Yj
d12 <- dXi1 + dpath[nbi1,nbj2] + d2Yj
d21 <- dXi2 + dpath[nbi2,nbj1] + d1Yj
d22 <- dXi2 + dpath[nbi2,nbj2] + d2Yj
d <- min(d11,d12,d21,d22)
}
# store result
crossdistmat[i,j] <- d
}
}
} else {
# C code
# convert indices to start at 0
from0 <- from - 1L
to0 <- to - 1L
Xsegmap <- Xpro - 1L
Ysegmap <- Ypro - 1L
DUP <- spatstat.options("dupC")
zz <- .C("lincrossdist",
np = as.integer(nX),
xp = as.double(P$x),
yp = as.double(P$y),
nq = as.integer(nY),
xq = as.double(Q$x),
yq = as.double(Q$y),
nv = as.integer(Lvert$n),
xv = as.double(Lvert$x),
yv = as.double(Lvert$y),
ns = as.double(L$n),
from = as.integer(from0),
to = as.integer(to0),
dpath = as.double(dpath),
psegmap = as.integer(Xsegmap),
qsegmap = as.integer(Ysegmap),
answer = as.double(crossdistmat),
DUP=DUP)
# PACKAGE="spatstat")
crossdistmat <- matrix(zz$answer, nX, nY)
}
return(crossdistmat)
}