https://github.com/cran/spatstat
Tip revision: 2910896e7be42bade87fb5b17a6239f86cb05760 authored by Adrian Baddeley on 15 June 2018, 06:48:03 UTC
version 1.56-0
version 1.56-0
Tip revision: 2910896
quadscheme.R
#
#
# quadscheme.S
#
# $Revision: 4.35 $ $Date: 2016/02/11 10:17:12 $
#
# quadscheme() generate a quadrature scheme from
# data and dummy point patterns.
#
# quadscheme.spatial() case where both patterns are unmarked
#
# quadscheme.replicated() case where data are multitype
#
#
#---------------------------------------------------------------------
quadscheme <- function(data, dummy, method="grid", ...) {
#
# generate a quadrature scheme from data and dummy patterns.
#
# Other arguments control how the quadrature weights are computed
#
data <- as.ppp(data)
if(missing(dummy)) {
# create dummy points
dummy <- default.dummy(data, method=method, ...)
# extract full set of parameters used to create dummy points
dp <- attr(dummy, "dummy.parameters")
# extract recommended parameters for computing weights
wp <- attr(dummy, "weight.parameters")
} else {
# user-supplied dummy points
if(!is.ppp(dummy)) {
# convert to ppp object
dummy <- as.ppp(dummy, data$window, check=FALSE)
# confine dummy points to data window
dummy <- dummy[data$window]
wp <- dp <- list()
} else {
# if it's already a ppp, it may have been created by default.dummy
dp <- attr(dummy, "dummy.parameters")
wp <- attr(dummy, "weight.parameters")
}
}
# arguments supplied directly to quadscheme()
# override any arguments passed as attributes
wp <- resolve.defaults(list(method=method), list(...), wp)
mX <- is.marked(data)
mD <- is.marked(dummy)
if(!mX && !mD)
Q <- do.call(quadscheme.spatial,
append(list(data, dummy, check=FALSE), wp))
else if(mX && !mD)
Q <- do.call(quadscheme.replicated,
append(list(data, dummy, check=FALSE), wp))
else if(!mX && mD)
stop("dummy points are marked but data are unmarked")
else
stop("marked data and marked dummy points -- sorry, this case is not implemented")
# record parameters used to make dummy points
Q$param$dummy <- dp
return(Q)
}
quadscheme.spatial <-
function(data, dummy, method=c("grid", "dirichlet"), ...) {
#
# generate a quadrature scheme from data and dummy patterns.
#
# The 'method' may be "grid" or "dirichlet"
#
# '...' are passed to gridweights() or dirichletWeights()
#
# quadscheme.spatial:
# for unmarked point patterns.
#
# weights are determined only by spatial locations
# (i.e. weight computations ignore any marks)
#
# No two points should have the same spatial location
#
check <- resolve.defaults(list(...), list(check=TRUE))$check
method <- match.arg(method)
data <- as.ppp(data, check=check)
dummy <- as.ppp(dummy, data$window, check=check)
# note data$window is the DEFAULT quadrature window
# applicable when 'dummy' does not contain a window
if(is.marked(data, dfok=TRUE))
warning("marks in data pattern - ignored")
if(is.marked(dummy, dfok=TRUE))
warning("marks in dummy pattern - ignored")
both <- as.ppp(concatxy(data, dummy), dummy$window, check=check)
switch(method,
grid={
w <- gridweights(both, window= dummy$window, ...)
},
dirichlet = {
w <- dirichletWeights(both, window=dummy$window, ...)
},
{
stop(paste("unrecognised method", sQuote(method)))
}
)
# parameters actually used to make weights
wp <- attr(w, "weight.parameters")
param <- list(weight = wp, dummy = NULL)
Q <- quad(data, dummy, w, param)
return(Q)
}
"quadscheme.replicated" <-
function(data, dummy, method=c("grid", "dirichlet"), ...) {
##
## generate a quadrature scheme from data and dummy patterns.
##
## The 'method' may be "grid" or "dirichlet"
##
## '...' are passed to gridweights() or dirichletWeights()
##
## quadscheme.replicated:
## for multitype point patterns.
##
## No two points in 'data'+'dummy' should have the same spatial location
check <- resolve.defaults(list(...), list(check=TRUE))$check
method <- match.arg(method)
data <- as.ppp(data, check=check)
dummy <- as.ppp(dummy, data$window, check=check)
## note data$window is the DEFAULT quadrature window
## unless otherwise specified in 'dummy'
ndata <- data$n
ndummy <- dummy$n
if(!is.marked(data))
stop("data pattern does not have marks")
if(is.marked(dummy, dfok=TRUE) && npoints(dummy) > 0)
warning("dummy points have marks --- ignored")
## first, ignore marks and compute spatial weights
P <- quadscheme.spatial(unmark(data), dummy, method, ...)
W <- w.quad(P)
iz <- is.data(P)
Wdat <- W[iz]
Wdum <- W[!iz]
## find the set of all possible marks
if(!is.multitype(data))
stop("data pattern is not multitype")
data.marks <- marks(data)
markset <- levels(data.marks)
nmarks <- length(markset)
## replicate dummy points, one copy for each possible mark
## -> dummy x {1,..,K}
dumdum <- cartesian(dummy, markset)
Wdumdum <- rep.int(Wdum, nmarks)
Idumdum <- rep.int(ndata + seq_len(ndummy), nmarks)
## also make dummy marked points at same locations as data points
## but with different marks
dumdat <- cartesian(unmark(data), markset)
Wdumdat <- rep.int(Wdat, nmarks)
Mdumdat <- marks(dumdat)
Idumdat <- rep.int(1:ndata, nmarks)
Mrepdat <- rep.int(data.marks, nmarks)
ok <- (Mdumdat != Mrepdat)
dumdat <- dumdat[ok,]
Wdumdat <- Wdumdat[ok]
Idumdat <- Idumdat[ok]
## combine the two dummy patterns
dumb <- superimpose(dumdum, dumdat, W=dummy$window, check=FALSE)
Wdumb <- c(Wdumdum, Wdumdat)
Idumb <- c(Idumdum, Idumdat)
## record the quadrature parameters
param <- list(weight = P$param$weight,
dummy = NULL,
sourceid=c(1:ndata, Idumb))
## wrap up
Q <- quad(data, dumb, c(Wdat, Wdumb), param)
return(Q)
}
"cartesian" <-
function(pp, markset, fac=TRUE) {
## given an unmarked point pattern 'pp'
## and a finite set of marks,
## create the marked point pattern which is
## the Cartesian product, consisting of all pairs (u,k)
## where u is a point of 'pp' and k is a mark in 'markset'
nmarks <- length(markset)
result <- ppp(rep.int(pp$x, nmarks),
rep.int(pp$y, nmarks),
window=pp$window,
check=FALSE)
marx <- rep.int(markset, rep.int(pp$n, nmarks))
if(fac)
marx <- factor(marx, levels=markset)
marks(result) <- marx
return(result)
}
validate.quad <- function(Q, fatal=FALSE, repair=TRUE, announce=FALSE) {
X <- Q$data
D <- Q$dummy
mX <- is.marked(X)
mD <- is.marked(D)
nbg <- function(whinge, fatal=FALSE, announce=FALSE) {
if(fatal)
stop(whinge, call.=FALSE)
else {
if(announce)
warning(whinge, call.=FALSE)
return(FALSE)
}
}
if(mX != mD) {
whinge <-
if(mX)
"data points are marked, but dummy points are not"
else
"dummy points are marked, but data points are not"
return(nbg(whinge, fatal, announce))
}
if(!mX)
return(TRUE)
# marked points
fX <- is.factor(Xmarx <- marks(X))
fD <- is.factor(Dmarx <- marks(D))
if(fX != fD) {
whinge <-
if(fX)
"data points are multitype, but dummy points are not"
else
"dummy points are multitype, but data points are not"
return(nbg(whinge, fatal, announce))
}
if(!fX)
return(TRUE)
# multitype points
lX <- levels(Xmarx)
lD <- levels(Dmarx)
if(length(lX) != length(lD) || any(lX != lD)) {
whinge <- "data and dummy points have different sets of possible marks"
return(nbg(whinge, fatal, announce))
}
return(TRUE)
}
pixelquad <- function(X, W=as.owin(X)) {
## make a quadscheme with a dummy point at every pixel
verifyclass(X, "ppp")
## convert window to mask if not already one
W <- as.owin(W)
M <- as.mask(W)
MM <- M$m
pixelarea <- M$xstep * M$ystep
## create pixel coordinates and corresponding row, column indices
rxy <- rasterxy.mask(M, drop=TRUE)
xx <- rxy$x
yy <- rxy$y
cc <- as.vector(col(MM)[MM])
rr <- as.vector(row(MM)[MM])
Nr <- M$dim[1]
Nc <- M$dim[2]
## dummy point pattern
dum <- ppp(xx, yy, window=W, check=FALSE)
## discretise data points
ij <- nearest.raster.point(X$x, X$y, M)
ijrow <- ij$row
ijcol <- ij$col
if(!is.marked(X)) {
## tabulate pixel locations of data points
Xtab <- table(row=factor(ijrow, levels=1:Nr),
col=factor(ijcol, levels=1:Nc))
## every pixel contains exactly one dummy point,
## so the total count of quadrature points in each pixel is:
Qtab <- Xtab + 1
## compute counting weights for data points
wdat <- 1/Qtab[cbind(ijrow, ijcol)]
## compute counting weights for dummy points
wdum <- 1/Qtab[cbind(rr, cc)]
} else {
marx <- marks(X)
## tabulate pixel locations and marks of data points
Xtab <- table(row=factor(ijrow, levels=1:Nr),
col=factor(ijcol, levels=1:Nc),
mark=marx)
## replicate dummy points (pixel centres) for each mark
dum <- cartesian(dum, levels(marx))
## every marked pixel contains exactly one dummy point,
## so the total count of quadrature points in each marked pixel is:
Qtab <- Xtab + 1
## compute counting weights for data points
wdat <- 1/Qtab[cbind(ijrow, ijcol, as.integer(marx))]
## compute counting weights for dummy points
nm <- length(levels(marx))
wdum <- 1/Qtab[cbind(rep.int(rr, nm),
rep.int(cc, nm),
rep(1:nm, each=length(rr)))]
}
## create quadrature scheme
wboth <- pixelarea * c(wdat, wdum)
Q <- quad(X, dum, wboth)
attr(Q, "M") <- M
return(Q)
}