Nothing
#' Test if case is orthogonal to segment.
#'
#' Diagnostic to check classification of case to a road segment.
#' @param case Numeric or Integer. Numeric ID of observed case.
#' @param segment Character. Segment ID. See \code{road.segments}.
#' @param observed Logical. \code{FALSE} observed case; \code{TRUE} simulated case (\code{regular.cases}).
#' @param coordinates Logical. Orthogonal projection coordinates.
#' @note This function is a diagnostic. It is not a guarantee of correct classification.
#' @return Logical \code{TRUE} or \code{FALSE}
#' @noRd
#' @examples
#' classifierAudit(case = 483, segment = "326-2")
#' plot(classifierAudit(case = 483, segment = "326-2"))
classifierAudit <- function(case = 483, segment = "326-2", observed = TRUE,
coordinates = FALSE) {
if (!is.numeric(case)) stop("case must be numeric.", call. = FALSE)
if (observed) {
if (case %in% unique(cholera::fatalities$case) == FALSE) {
stop("Observed case must be a whole number between 1 and 578.",
call. = FALSE)
}
obs <- cholera::fatalities[cholera::fatalities$case == case, c("x", "y")]
} else {
if (case %in% seq_len(nrow(cholera::regular.cases)) == FALSE) {
stop("Simulated case must be a whole number between 1 and 4993.",
call. = FALSE)
}
obs <- cholera::regular.cases[case, c("x", "y")]
}
seg.data <- cholera::road.segments[cholera::road.segments$id == segment,
c("x1", "y1", "x2", "y2")]
seg.df <- data.frame(x = c(seg.data$x1, seg.data$x2),
y = c(seg.data$y1, seg.data$y2))
ols <- stats::lm(y ~ x, data = seg.df)
segment.slope <- stats::coef(ols)[2]
segment.intercept <- stats::coef(ols)[1]
orthogonal.slope <- -1 / segment.slope
orthogonal.intercept <- obs$y - orthogonal.slope * obs$x
x.proj <- (orthogonal.intercept - segment.intercept) /
(segment.slope - orthogonal.slope)
y.proj <- segment.slope * x.proj + segment.intercept
if (coordinates) {
data.frame(x.proj, y.proj, row.names = NULL)
} else {
# Bisection / Intersection test
distB <- stats::dist(rbind(seg.df[1, ], c(x.proj, y.proj))) +
stats::dist(rbind(seg.df[2, ], c(x.proj, y.proj)))
distance <- stats::dist(seg.df)
proj <- rbind(obs, data.frame(x = x.proj, y = y.proj, row.names = NULL))
ortho.distance <- stats::dist(proj)
out <- list(case = case,
segment = segment,
ols = ols,
seg.df = seg.df,
obs = obs,
distance = ortho.distance,
test = signif(distance) == signif(distB),
coords = data.frame(x.proj, y.proj, row.names = NULL))
class(out) <- "classifier_audit"
out
}
}
#' Return result of classifierAudit().
#'
#' @param x An object of class "classifier_audit" created by \code{classifierAudit()}.
#' @param ... Additional parameters.
#' @return An R data frame.
#' @noRd
#' @examples
#' classifierAudit(case = 483, segment = "326-2")
#' print(classifierAudit(case = 483, segment = "326-2"))
print.classifier_audit <- function(x, ...) {
if (!inherits(x, "classifier_audit")) {
stop('"x"\'s class needs to be "classifier_audit".')
}
print(x$test)
}
#' Plot result of classifierAudit().
#'
#' Plot case, segment and orthogonal projector.
#' @param x An object of class "classifier_audit" created by \code{classifierAudit()}.
#' @param zoom Logical or Numeric. A numeric value >= 0 controls the degree of zoom. The default is 0.5.
#' @param unit Character. Unit of distance: "meter" (the default), "yard" or "native". "native" returns the map's native scale. "unit" is meaningful only when "weighted" is \code{TRUE}. See \code{vignette("roads")} for information on unit distances.
#' @param ... Additional parameters.
#' @return A base R graphic.
#' @noRd
#' @examples
#' plot(classifierAudit(case = 483, segment = "326-2"))
plot.classifier_audit <- function(x, zoom = 0.5, unit = "meter", ...) {
if (unit %in% c("meter", "yard", "native") == FALSE) {
stop('unit must be "meter", "yard" or "native".')
}
obs <- x$obs
coords <- x$coords
x.proj <- coords$x.proj
y.proj <- coords$y.proj
segmentLocator(x$segment, zoom = zoom, add.title = FALSE,
add.subtitle = FALSE)
# Bisection / Intersection test
distB <- stats::dist(rbind(x$seg.df[1, ], c(x.proj, y.proj))) +
stats::dist(rbind(x$seg.df[2, ], c(x.proj, y.proj)))
if (x$test) {
points(obs, col = "green")
points(x.proj, y.proj, pch = 0, col = "green")
arrows(obs$x, obs$y, x.proj, y.proj, length = 0.1, col = "green", lwd = 3)
} else {
points(obs, col = "dodgerblue")
points(x.proj, y.proj, pch = 0, col = "dodgerblue")
arrows(obs$x, obs$y, x.proj, y.proj, length = 0.1, col = "dodgerblue",
lwd = 3)
}
if (unit == "meter") {
ortho.dist <- unitMeter(x$distance, "meter")
d.unit <- "m"
} else if (unit == "yard") {
ortho.dist <- unitMeter(x$distance, "yard")
d.unit <- "yd"
} else if (unit == "native") {
ortho.dist <- unitMeter(x$distance, "native")
d.unit <- "units"
}
nm <- cholera::road.segments[cholera::road.segments$id == x$segment, "name"]
title(main = paste0(nm, ": Segment # ", x$segment, ", Case # ", x$case),
sub = paste(x$test, "; ortho.dist =", round(ortho.dist, 1), d.unit))
}
#' Visually check orthogonal projection of case fatalities to road segments (address).
#'
#' @param case Numeric. case ID from \code{fatalities}.
#' @param rd.seg Character. Road segment ID.
#' @noRd
orthoProjector <- function(case = 483, rd.seg = "326-2") {
case.data <- cholera::fatalities[cholera::fatalities$case == case, ]
seg.data <- cholera::road.segments[cholera::road.segments$id == rd.seg,
c("x1", "y1", "x2", "y2")]
seg.df <- data.frame(x = c(seg.data$x1, seg.data$x2),
y = c(seg.data$y1, seg.data$y2))
ols <- stats::lm(y ~ x, data = seg.df)
segment.slope <- stats::coef(ols)[2]
segment.intercept <- stats::coef(ols)[1]
orthogonal.slope <- -1 / segment.slope
orthogonal.intercept <- case.data$y - orthogonal.slope * case.data$x
x.proj <- (orthogonal.intercept - segment.intercept) /
(segment.slope - orthogonal.slope)
y.proj <- segment.slope * x.proj + segment.intercept
x.rng <- range(seg.data[, c("x1", "x2")], case.data$x, x.proj)
y.rng <- range(seg.data[, c("y1", "y2")], case.data$y, y.proj)
distB <- stats::dist(rbind(seg.df[1, ], c(x.proj, y.proj))) +
stats::dist(rbind(seg.df[2, ], c(x.proj, y.proj)))
bisect.test <- signif(stats::dist(seg.df)) == signif(distB)
plot(seg.data[, c("x1", "y1")], xlim = x.rng, ylim = y.rng, xlab = "x",
ylab = "y", asp = 1)
points(seg.data[, c("x2", "y2")])
segments(seg.data$x1, seg.data$y1, seg.data$x2, seg.data$y2)
text(case.data$x, case.data$y, pos = 1, labels = case)
title(main = rd.seg)
if (bisect.test) {
points(case.data$x, case.data$y, col = "green", pch = 16)
arrows(case.data$x, case.data$y, x.proj, y.proj, length = 0.1,
col = "green")
} else {
points(case.data$x, case.data$y, col = "red", pch = 16)
segments(seg.data$x2, seg.data$y2, x.proj, y.proj, lty = "dotted",
col = "gray")
arrows(case.data$x, case.data$y, x.proj, y.proj, length = 0.1, col = "red")
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.