#' Compute network graph of roads, cases and pumps.
#'
#' Assembles cases, pumps and road into a network graph.
#' @param vestry Logical. Use Vestry Report pump data.
#' @param case.set Character. "observed", "expected", or "snow". "snow" captures John Snow's annotation of the Broad Street pump neighborhood printed in the Vestry report version of the map.
#' @param embed Logical. Embed cases and pumps into road network.
#' @param embed.landmarks Logical. Embed landmarks into road network.
#' @export
#' @return An R list of nodes, edges and an 'igraph' network graph.
neighborhoodData <- function(vestry = FALSE, case.set = "observed",
embed = TRUE, embed.landmarks = TRUE) {
if (case.set %in% c("observed", "expected", "snow") == FALSE) {
stop('case.set must be "observed", "expected" or "snow".')
}
args <- list(embed = embed, embed.landmarks = embed.landmarks,
vestry = vestry, case.set = case.set)
node.data <- do.call("nodeData", args)
nodes <- node.data$nodes
edges <- node.data$edges
g <- node.data$g
nodes.pump <- nodes[nodes$pump != 0, ]
nodes.pump <- nodes.pump[order(nodes.pump$pump), c("pump", "node")]
out <- list(g = g, nodes = nodes, edges = edges, nodes.pump = nodes.pump)
class(out) <- "neighborhood_data"
out
}
#' Plot method for neighborhoodData().
#'
#' Visualize underlying road network (with or without cases and pumps).
#' @param x An 'igraph' object of class "neighborhood_data" created by \code{neighborhoodData()}.
#' @param ... Additional plotting parameters.
#' @return A base R plot.
#' @export
#' @examples
#' \donttest{
#' plot(neighborhoodData())
#' plot(neighborhoodData(embed = FALSE))
#' }
plot.neighborhood_data <- function(x, ...) {
plot(x$g, vertex.label = NA, vertex.size = 2, ...)
}
## auxilliary functions ##
nodeData <- function(embed = TRUE, embed.landmarks = FALSE, vestry = FALSE,
case.set = "observed") {
if (embed) {
if (vestry) {
ortho.pump <- cholera::ortho.proj.pump.vestry
} else {
ortho.pump <- cholera::ortho.proj.pump
}
if (case.set == "observed") {
sel <- cholera::ortho.proj$case %in% cholera::fatalities.address$anchor
case.segments <- unique(cholera::ortho.proj[sel, "road.segment"])
} else if (case.set == "expected") {
sim.proj <- cholera::sim.ortho.proj
case.segments <- unique(sim.proj$road.segment)
case.segments <- case.segments[is.na(case.segments) == FALSE]
} else if (case.set == "snow") {
sel <- cholera::ortho.proj$case %in% cholera::snow.neighborhood
case.segments <- unique(cholera::ortho.proj[sel, "road.segment"])
}
if (embed.landmarks) {
landmark.segments <- cholera::landmarks$road.segment
site.segments <- unique(c(case.segments, landmark.segments))
} else {
site.segments <- case.segments
}
site.pump <- intersect(ortho.pump$road.segment, site.segments)
site.no_pump <- setdiff(site.segments, ortho.pump$road.segment)
no_site.pump <- setdiff(ortho.pump$road.segment, site.segments)
edits <- c(site.pump, site.no_pump, no_site.pump)
if (case.set %in% c("observed", "snow")) {
if (vestry) {
nodes <- lapply(edits, embedSites, ortho.pump, vestry = TRUE)
edges <- lapply(edits, embedSites, ortho.pump, type = "edges",
vestry = TRUE)
} else {
nodes <- lapply(edits, embedSites, ortho.pump)
edges <- lapply(edits, embedSites, ortho.pump, type = "edges")
}
} else if (case.set == "expected") {
if (vestry) {
nodes <- lapply(edits, embedSites, ortho.pump, case.set = "expected",
vestry = TRUE)
edges <- lapply(edits, embedSites, ortho.pump, type = "edges",
case.set = "expected", vestry = TRUE)
} else {
nodes <- lapply(edits, embedSites, ortho.pump, case.set = "expected")
edges <- lapply(edits, embedSites, ortho.pump, type = "edges",
case.set = "expected")
}
}
# Edges #
edges <- do.call(rbind, edges)
road.segments <- cholera::road.segments
road.segments <- road.segments[road.segments$id %in% edits == FALSE, ]
road.segments$node1 <- paste0(road.segments$x1, "_&_", road.segments$y1)
road.segments$node2 <- paste0(road.segments$x2, "_&_", road.segments$y2)
if (case.set %in% c("observed", "snow")) {
road.segments$id2 <- paste0(road.segments$id, "a")
} else if (case.set == "expected") {
road.segments$id2 <- paste0(road.segments$id, "-1")
}
edges <- rbind(edges, road.segments)
edges$d <- sqrt((edges$x1 - edges$x2)^2 + (edges$y1 - edges$y2)^2)
# Nodes #
nodes <- do.call(rbind, nodes)
select1 <- road.segments$node1 %in% nodes$node == FALSE
select2 <- road.segments$node2 %in% nodes$node == FALSE
endpt.node1 <- road.segments[select1, "node1"]
endpt.node2 <- road.segments[select2, "node2"]
null.nodes <- lapply(union(endpt.node1, endpt.node2), function(x) {
sel <- road.segments$node1 == x | road.segments$node2 == x
dat <- road.segments[sel, ]
endpt1 <- dat$node1 == x & dat$node2 != x
endpt2 <- dat$node1 != x & dat$node2 == x
if (any(endpt1)) {
data.sel <- unique(dat[endpt1, c("x1", "y1")])
} else if (any(endpt2)) {
data.sel <- unique(dat[endpt2, c("x2", "y2")])
}
data.sel <- stats::setNames(data.sel, c("x.proj", "y.proj"))
data.frame(data.sel, anchor = 0, pump = 0, node = x)
})
null.nodes <- do.call(rbind, null.nodes)
nodes <- rbind(nodes, null.nodes)
# Network Graph #
edge.list <- edges[, c("node1", "node2")]
g <- igraph::graph_from_data_frame(edge.list, directed = FALSE)
list(nodes = nodes, edges = edges, g = g)
} else {
road.segments$node1 <- paste0(road.segments$x1, "_&_", road.segments$y1)
road.segments$node2 <- paste0(road.segments$x2, "_&_", road.segments$y2)
if (case.set == "observed") {
road.segments$id2 <- paste0(road.segments$id, "a")
} else if (case.set == "expected") {
road.segments$id2 <- paste0(road.segments$id, "-1")
}
road.segments$d <- sqrt((road.segments$x1 - road.segments$x2)^2 +
(road.segments$y1 - road.segments$y2)^2)
rs1 <- road.segments[, c("x1", "y1", "node1")]
rs2 <- road.segments[, c("x2", "y2", "node2")]
a <- stats::setNames(rs1, c("x.proj", "y.proj", "node"))
b <- stats::setNames(rs2, c("x.proj", "y.proj", "node"))
nodes <- unique(rbind(a, b))
nodes <- data.frame(nodes[, c("x.proj", "y.proj")],
anchor = 0,
pump = 0,
node = nodes$node,
stringsAsFactors = FALSE)
edge.list <- road.segments[, c("node1", "node2")]
g <- igraph::graph_from_data_frame(edge.list, directed = FALSE)
list(nodes = nodes, edges = road.segments, g = g)
}
}
embedSites <- function(id, ortho.pump, type = "nodes", case.set = "observed",
vestry = FALSE) {
road.data <- cholera::road.segments[cholera::road.segments$id == id, ]
## embed filters ##
CaseObs <- id %in% cholera::ortho.proj$road.segment
Landmark <- id %in% cholera::landmarks$road.segment
Pump <- id %in% ortho.pump$road.segment
if (case.set == "expected") {
sim.proj <- cholera::sim.ortho.proj
CaseExp <- id %in% sim.proj$road.segment
}
if (case.set == "observed") {
if (CaseObs) {
sel <- cholera::ortho.proj$road.segment %in% id
road.fatalities <- cholera::ortho.proj[sel, ]
sel <- road.fatalities$case %in% cholera::fatalities.address$anchor
cases <- road.fatalities[sel, "case"]
road.address <- road.fatalities[road.fatalities$case %in% cases, ]
case.data <- data.frame(road.address[, c("x.proj", "y.proj")],
anchor = road.address$case,
pump = 0)
}
} else if (case.set == "expected") {
if (CaseExp) {
road.fatalities <- sim.proj[sim.proj$road.segment %in% id, ]
sel <- road.fatalities$case[road.fatalities$case %in% sim.proj$case]
road.address <- road.fatalities[road.fatalities$case %in% sel, ]
case.data <- data.frame(road.address[, c("x.proj", "y.proj")],
anchor = road.address$case,
pump = 0)
}
}
if (Landmark) {
sel <- cholera::landmarks$road.segment %in% id
road.landmarks <- cholera::landmarks[sel, ]
landmark.data <- data.frame(road.landmarks[, c("x.proj", "y.proj")],
anchor = road.landmarks$case,
pump = 0)
}
endptA <- data.frame(x.proj = road.data$x1,
y.proj = road.data$y1,
anchor = 0,
pump = 0)
endptB <- data.frame(x.proj = road.data$x2,
y.proj = road.data$y2,
anchor = 0,
pump = 0)
if (case.set == "observed") {
if (CaseObs & Landmark) {
site.seg <- CaseObs & Landmark
rds <- rbind(case.data, landmark.data)
} else if (CaseObs & !Landmark) {
site.seg <- CaseObs
rds <- case.data
} else if (!CaseObs & Landmark) {
site.seg <- Landmark
rds <- landmark.data
} else site.seg <- FALSE
} else if (case.set == "expected") {
if (CaseExp & Landmark) {
site.seg <- CaseExp & Landmark
rds <- rbind(case.data, landmark.data)
} else if (CaseExp & !Landmark) {
site.seg <- CaseExp
rds <- case.data
} else if (!CaseExp & Landmark) {
site.seg <- Landmark
rds <- landmark.data
} else site.seg <- FALSE
}
if (Pump) {
pump.data <- ortho.pump[ortho.pump$road.segment == id, ]
ps <- data.frame(pump.data[, c("x.proj", "y.proj")],
anchor = 0,
pump = pump.data$pump.id)
}
if (site.seg & Pump) {
dat <- rbind(rds, ps)
nodes <- rbind(endptA, dat, endptB)
nodes <- nodes[order(nodes$x.proj), ]
} else if (site.seg & !Pump) {
if (case.set == "observed") {
pts <- rbind(endptA, rds, endptB)
duplicates <- duplicateNode(pts) # site at intersection
if (any(duplicates)) {
if (which(duplicates) == 1) {
nodes <- rbind(rds, endptB)
} else if (which(duplicates) == 2) {
nodes <- rbind(endptA, rds)
}
} else {
nodes <- rbind(endptA, rds, endptB)
}
} else if (case.set == "expected") {
nodes <- rbind(endptA, rds, endptB)
}
nodes <- nodes[order(nodes$x.proj), ]
} else if (!site.seg & Pump) {
nodes <- rbind(endptA, ps, endptB)
nodes <- nodes[order(nodes$x.proj), ]
} else {
nodes <- rbind(endptA, endptB)
nodes <- nodes[order(nodes$x.proj), ]
}
row.names(nodes) <- NULL
nodes$node <- paste0(nodes$x.proj, "_&_", nodes$y.proj)
# Edges
sel <- c("x.proj", "y.proj")
edges <- cbind(nodes[-nrow(nodes), sel], nodes[-1, sel])
names(edges) <- c("x1", "y1", "x2", "y2")
edges$node1 <- paste0(edges$x1, "_&_", edges$y1)
edges$node2 <- paste0(edges$x2, "_&_", edges$y2)
edges <- cbind(road.data[1, c("street", "id", "name")], edges,
row.names = NULL)
if (case.set == "observed") {
edges$id2 <- paste0(edges$id, letters[seq_len(nrow(edges))])
} else if (case.set == "expected") {
edges$id2 <- paste0(edges$id, "-", seq_len(nrow(edges)))
}
if (type == "nodes") {
nodes
} else if (type == "edges") {
edges
}
}
duplicateNode <- function(dat, sig.fig = 5L) {
vars <- c("x.proj", "y.proj")
idx <- seq_len(nrow(dat[-1, ]))
idx <- list(idx, idx + 1)
vapply(idx, function(sel) {
all(round(dat[sel[1], vars], sig.fig) == round(dat[sel[2], vars], sig.fig))
}, logical(1L))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.