#####################
# Function chooseCN
#####################
#' Function to choose a connection network
#'
#' The function \code{chooseCN} is a simple interface to build a connection
#' network (CN) from xy coordinates. The user chooses from 6 types of graph and
#' one additional weighting scheme. \code{chooseCN} calls functions from
#' appropriate packages, handles non-unique coordinates and returns a
#' connection network either with classe \code{nb} or \code{listw}. For graph
#' types 1-4, duplicated locations are not accepted and will issue an error.
#'
#' There are 7 kinds of graphs proposed: \cr Delaunay triangulation (type 1)\cr
#' Gabriel graph (type 2)\cr Relative neighbours (type 3)\cr Minimum spanning
#' tree (type 4)\cr Neighbourhood by distance (type 5)\cr K nearests neighbours
#' (type 6)\cr Inverse distances (type 7)\cr
#'
#' The last option (type=7) is not a true neighbouring graph: all sites are
#' neighbours, but the spatial weights are directly proportional to the
#' inversed spatial distances.\cr Also not that in this case, the output of the
#' function is always a \code{listw} object, even if \code{nb} was
#' requested.\cr
#'
#' The choice of the connection network has been discuted on the adegenet
#' forum. Please search the archives from adegenet website (section 'contact')
#' using 'graph' as keyword.
#'
#' @param xy an matrix or data.frame with two columns for x and y coordinates.
#' @param ask a logical stating whether graph should be chosen interactively
#' (TRUE,default) or not (FALSE). Set to FALSE if \code{type} is provided.
#' @param type an integer giving the type of graph (see details).
#' @param result.type a character giving the class of the returned object.
#' Either "nb" (default) or "listw", both from \code{spdep} package. See
#' details.
#' @param d1 the minimum distance between any two neighbours. Used if
#' \code{type=5.}
#' @param d2 the maximum distance between any two neighbours. Used if
#' \code{type=5}. Can also be a character: "dmin" for the minimum distance so
#' that each site has at least one connection, or "dmax" to have all sites
#' connected (despite the later has no sense).
#' @param k the number of neighbours per point. Used if \code{type=6}.
#' @param a the exponent of the inverse distance matrix. Used if \code{type=7}.
#' @param dmin the minimum distance between any two distinct points. Used to
#' avoid infinite spatial proximities (defined as the inversed spatial
#' distances). Used if \code{type=7}.
#' @param plot.nb a logical stating whether the resulting graph should be
#' plotted (TRUE, default) or not (FALSE).
#' @param edit.nb a logical stating whether the resulting graph should be
#' edited manually for corrections (TRUE) or not (FALSE, default).
#' @param check.duplicates a logical indicating if duplicate coordinates should be detected; this can be an issue for some graphs; TRUE by default.
#'
#' @return Returns a connection network having the class \code{nb} or
#' \code{listw}. The xy coordinates are passed as attribute to the created
#' object.
#'
#' @author Thibaut Jombart \email{t.jombart@@imperial.ac.uk}
#'
#' @seealso \code{\link{spca}}
#' @keywords spatial utilities
#' @examples
#'
#' \dontrun{
#' data(nancycats)
#'
#' par(mfrow=c(2,2))
#' cn1 <- chooseCN(nancycats@@other$xy,ask=FALSE,type=1)
#' cn2 <- chooseCN(nancycats@@other$xy,ask=FALSE,type=2)
#' cn3 <- chooseCN(nancycats@@other$xy,ask=FALSE,type=3)
#' cn4 <- chooseCN(nancycats@@other$xy,ask=FALSE,type=4)
#' par(mfrow=c(1,1))
#' }
#'
#' @export chooseCN
#' @import ade4
# importFrom spdep "listw2mat" "lag.listw" "card"
#'
chooseCN <- function(xy, ask = TRUE, type = NULL, result.type = "nb",
d1 = NULL, d2 = NULL, k = NULL, a = NULL,
dmin = NULL, plot.nb = TRUE, edit.nb = FALSE,
check.duplicates = TRUE){
if (!requireNamespace("spdep", quietly=TRUE)) {
install <- paste0('install.packages(', shQuote("spdep"), ")")
msg <- c("The spdep package is required. Please use `", install, "` to install it")
stop(paste(msg, collapse = ""))
}
if(is.data.frame(xy)) xy <- as.matrix(xy)
if(ncol(xy) != 2) stop("xy does not have two columns.")
if(any(is.na(xy))) stop("NA entries in xy.")
result.type <- tolower(result.type)
if(is.null(type) & !ask) stop("Non-interactive mode but no graph chosen; please provide a value for 'type' argument.")
res <- list()
if(!is.null(d2)){
if(d2=="dmin"){
tempmat <- as.matrix(dist(xy))
d2min <- max(apply(tempmat, 1, function(r) min(r[r>1e-12])))
d2min <- d2min * 1.0001 # to avoid exact number problem
d2 <- d2min
} else if(d2=="dmax"){
d2max <- max(dist(xy))
d2max <- d2max * 1.0001 # to avoid exact number problem
d2 <- d2max
}
} # end handle d2
d1.first <- d1
d2.first <- d2
k.first <- k
## handle type argument
if(!is.null(type)){
type <- as.integer(type)
if(type < 1 |type > 7) stop("type must be between 1 and 7")
ask <- FALSE
}
## check for uniqueness of coordinates
if(check.duplicates && any(xyTable(xy)$number>1)){ # if duplicate coords
DUPLICATE.XY <- TRUE
} else {
DUPLICATE.XY <- FALSE
}
## if(is.null(type) & !ask) { type <- 1 }
### begin large while ###
chooseAgain <- TRUE
while(chooseAgain){
# re-initialisation of some variables
d1 <- d1.first
d2 <- d2.first
k <- k.first
## read type from console
if(ask){
temp <- TRUE
while(temp){
cat("\nChoose a connection network:\n")
cat("\t Delaunay triangulation (type 1)\n")
cat("\t Gabriel graph (type 2)\n")
cat("\t Relative neighbours (type 3)\n")
cat("\t Minimum spanning tree (type 4)\n")
cat("\t Neighbourhood by distance (type 5)\n")
cat("\t K nearest neighbours (type 6)\n")
cat("\t Inverse distances (type 7)\n")
cat("Answer: ")
type <- as.integer(readLines(con = getOption('adegenet.testcon'), n = 1))
temp <- type < 1 |type > 7
if(temp) cat("\nWrong answer\n")
if(type %in% 1:4 & DUPLICATE.XY){
cat("\n\n== PROBLEM DETECTED ==")
cat("\nDuplicate locations detected\nPlease choose another graph (5-7) or add random noise to locations (see ?jitter).\n")
temp <- TRUE
}
} # end while
}
##
## warning about duplicate xy coords
if(type %in% 1:4 & DUPLICATE.XY){
stop("Duplicate locations detected and incompatible with graph type 1-4.\nPlease choose another graph (5-7) or add random noise to locations (see ?jitter).")
}
## graph types
## type 1: Delaunay
if(type==1){
cn <- spdep::tri2nb(xy)
}
# type 2: Gabriel
if(type==2){
cn <- spdep::gabrielneigh(xy)
cn <- spdep::graph2nb(cn, sym=TRUE)
}
## type 3: Relative neighbours
if(type==3){
cn <- spdep::relativeneigh(xy)
cn <- spdep::graph2nb(cn, sym=TRUE)
}
## type 4: Minimum spanning tree
if(type==4){
cn <- ade4::mstree(dist(xy)) # there is also a spdep::mstree
cn <- neig2nb(cn)
}
## type 5: Neighbourhood by distance
if(type==5){
if(is.null(d1) |is.null(d2)){
tempmat <- as.matrix(dist(xy))
d2min <- max(apply(tempmat, 1, function(r) min(r[r>1e-12])))
d2min <- d2min * 1.0001 # to avoid exact number problem
d2max <- max(dist(xy))
d2max <- d2max * 1.0001 # to avoid exact number problem
dig <- options("digits")
options("digits=5")
cat("\n Enter minimum distance: ")
d1 <- as.numeric(readLines(con = getOption('adegenet.testcon'), n = 1))
cat("\n Enter maximum distance \n(dmin=", d2min, ", dmax=", d2max, "): ")
d2 <- readLines(con = getOption('adegenet.testcon'), n = 1)
## handle character
if(d2=="dmin") {
d2 <- d2min
} else if(d2=="dmax") {
d2 <- d2max
} else {
d2 <- as.numeric(d2)
}
## restore initial digit option
options(dig)
}
# avoid that a point is its neighbour
dmin <- mean(dist(xy))/100000
if(d1<dmin) d1 <- dmin
if(d2<d1) stop("d2 < d1")
cn <- spdep::dnearneigh(x=xy, d1=d1, d2=d2)
}
## type 6: K nearests
if(type==6){
if(is.null(k)) {
cat("\n Enter the number of neighbours: ")
k <- as.numeric(readLines(con = getOption('adegenet.testcon'), n = 1))
}
cn <- spdep::knearneigh(x=xy, k=k)
cn <- spdep::knn2nb(cn, sym=TRUE)
}
## type 7: inverse distances
if(type==7){
if(is.null(a)) {
cat("\n Enter the exponent: ")
a <- as.numeric(readLines(con = getOption('adegenet.testcon'), n = 1))
}
cn <- as.matrix(dist(xy))
if(is.null(dmin)) {
cat("\n Enter the minimum distance \n(range = 0 -", max(cn),"): ")
dmin <- as.numeric(readLines(con = getOption('adegenet.testcon'), n = 1))
}
if(a<1) { a <- 1 }
thres <- mean(cn)/1e8
if(dmin < thres) dmin <- thres
cn[cn < dmin] <- dmin
cn <- 1/(cn^a)
diag(cn) <- 0
cn <- prop.table(cn,1)
plot.nb <- FALSE
edit.nb <- FALSE
result.type <- "listw"
} # end type 7
## end graph types
if(ask & plot.nb) {
plot(cn,xy)
cat("\nKeep this graph (y/n)? ")
ans <- tolower(readLines(con = getOption('adegenet.testcon'), n=1))
if(ans=="n") {chooseAgain <- TRUE} else {chooseAgain <- FALSE}
}
else if(plot.nb){
plot(cn,xy)
chooseAgain <- FALSE
}
else {chooseAgain <- FALSE}
}
### end large while
if(edit.nb) {cn <- edit(cn,xy)}
if(result.type == "listw") {
if(type!=7) {
cn <- spdep::nb2listw(cn, style="W", zero.policy=TRUE)
} else {
cn <- spdep::mat2listw(cn)
cn$style <- "W"
}
}
res <- cn
attr(res,"xy") <- xy
return(res)
} # end chooseCN
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.