diff --git a/.Rbuildignore b/.Rbuildignore
index e4bff0b18..f1a8c3636 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -16,6 +16,7 @@
 ^data-raw$
 ^vignettes/articles$
 ^.git-blame-ignore-revs$
+^DEVELOPMENT\.md$
 ^doc$
 ^Meta$
 ^.lintr$
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index d21127830..8dbe7ea7f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,7 +23,7 @@ URL: https://github.com/cmu-delphi/epipredict/,
     https://cmu-delphi.github.io/epipredict
 BugReports: https://github.com/cmu-delphi/epipredict/issues/
 Depends:
-    epiprocess (>= 0.6.0),
+    epiprocess (>= 0.7.5),
     parsnip (>= 1.0.0),
     R (>= 3.5.0)
 Imports:
@@ -32,6 +32,7 @@ Imports:
     distributional,
     dplyr,
     generics,
+    ggplot2,
     glue,
     hardhat (>= 1.3.0),
     magrittr,
@@ -51,7 +52,6 @@ Suggests:
     data.table,
     epidatr (>= 1.0.0),
     fs,
-    ggplot2,
     knitr,
     lubridate,
     poissonreg,
diff --git a/NAMESPACE b/NAMESPACE
index 969daf0c0..bfc4b696d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,6 +10,8 @@ S3method(adjust_frosting,frosting)
 S3method(apply_frosting,default)
 S3method(apply_frosting,epi_workflow)
 S3method(augment,epi_workflow)
+S3method(autoplot,canned_epipred)
+S3method(autoplot,epi_workflow)
 S3method(bake,check_enough_train_data)
 S3method(bake,epi_recipe)
 S3method(bake,step_epi_ahead)
@@ -23,6 +25,7 @@ S3method(detect_layer,workflow)
 S3method(epi_keys,data.frame)
 S3method(epi_keys,default)
 S3method(epi_keys,epi_df)
+S3method(epi_keys,epi_workflow)
 S3method(epi_keys,recipe)
 S3method(epi_recipe,default)
 S3method(epi_recipe,epi_df)
@@ -128,6 +131,7 @@ export(arx_class_epi_workflow)
 export(arx_classifier)
 export(arx_fcast_epi_workflow)
 export(arx_forecaster)
+export(autoplot)
 export(bake)
 export(cdc_baseline_args_list)
 export(cdc_baseline_forecaster)
@@ -215,6 +219,7 @@ importFrom(dplyr,ungroup)
 importFrom(epiprocess,growth_rate)
 importFrom(generics,augment)
 importFrom(generics,fit)
+importFrom(ggplot2,autoplot)
 importFrom(hardhat,refresh_blueprint)
 importFrom(hardhat,run_mold)
 importFrom(magrittr,"%>%")
diff --git a/R/autoplot.R b/R/autoplot.R
new file mode 100644
index 000000000..fa3f26269
--- /dev/null
+++ b/R/autoplot.R
@@ -0,0 +1,298 @@
+#' @importFrom ggplot2 autoplot
+#' @export
+ggplot2::autoplot
+
+#' Automatically plot an `epi_workflow` or `canned_epipred` object
+#'
+#' For a fit workflow, the training data will be displayed, the response by
+#' default. If `predictions` is not `NULL` then point and interval forecasts
+#' will be shown as well. Unfit workflows will result in an error, (you
+#' can simply call `autoplot()` on the original `epi_df`).
+#'
+#'
+#'
+#'
+#' @inheritParams epiprocess::autoplot.epi_df
+#' @param object An `epi_workflow`
+#'
+#' @param predictions A data frame with predictions. If `NULL`, only the
+#'   original data is shown.
+#' @param .levels A numeric vector of levels to plot for any prediction bands.
+#'   More than 3 levels begins to be difficult to see.
+#' @param ... Ignored
+#' @param .color_by A character string indicating how to color the data. See
+#'   `epiprocess::autoplot.epi_df()` for more details.
+#' @param .facet_by A character string indicating how to facet the data. See
+#'  `epiprocess::autoplot.epi_df()` for more details.
+#' @param .base_color If available, prediction bands will be shown with this
+#'   color.
+#' @param .point_pred_color If available, point forecasts will be shown with this
+#'   color.
+#' @param .max_facets The maximum number of facets to show. If the number of
+#'  facets is greater than this value, only the top facets will be shown.
+#'
+#' @name autoplot-epipred
+#' @examples
+#' jhu <- case_death_rate_subset %>%
+#'   filter(time_value >= as.Date("2021-11-01"))
+#'
+#' r <- epi_recipe(jhu) %>%
+#'   step_epi_lag(death_rate, lag = c(0, 7, 14)) %>%
+#'   step_epi_ahead(death_rate, ahead = 7) %>%
+#'   step_epi_lag(case_rate, lag = c(0, 7, 14)) %>%
+#'   step_epi_naomit()
+#'
+#' f <- frosting() %>%
+#'   layer_residual_quantiles(
+#'     quantile_levels = c(.025, .1, .25, .75, .9, .975)
+#'   ) %>%
+#'   layer_threshold(dplyr::starts_with(".pred")) %>%
+#'   layer_add_target_date()
+#'
+#' wf <- epi_workflow(r, parsnip::linear_reg(), f) %>% fit(jhu)
+#'
+#' autoplot(wf)
+#'
+#' latest <- jhu %>% dplyr::filter(time_value >= max(time_value) - 14)
+#' preds <- predict(wf, latest)
+#' autoplot(wf, preds, .max_facets = 4)
+#'
+#' # ------- Show multiple horizons
+#'
+#' p <- lapply(c(7, 14, 21, 28), \(h) {
+#'   r <- epi_recipe(jhu) %>%
+#'     step_epi_lag(death_rate, lag = c(0, 7, 14)) %>%
+#'     step_epi_ahead(death_rate, ahead = h) %>%
+#'     step_epi_lag(case_rate, lag = c(0, 7, 14)) %>%
+#'     step_epi_naomit()
+#'   ewf <- epi_workflow(r, parsnip::linear_reg(), f) %>% fit(jhu)
+#'   td <- get_test_data(r, jhu)
+#'   predict(ewf, new_data = td)
+#' })
+#'
+#' p <- do.call(rbind, p)
+#' autoplot(wf, p, .max_facets = 4)
+#'
+#' # ------- Plotting canned forecaster output
+#'
+#' jhu <- case_death_rate_subset %>% filter(time_value >= as.Date("2021-11-01"))
+#' flat <- flatline_forecaster(jhu, "death_rate")
+#' autoplot(flat, .max_facets = 4)
+#'
+#' arx <- arx_forecaster(jhu, "death_rate", c("case_rate", "death_rate"),
+#'   args_list = arx_args_list(ahead = 14L)
+#' )
+#' autoplot(arx, .max_facets = 6)
+NULL
+
+#' @export
+#' @rdname autoplot-epipred
+autoplot.epi_workflow <- function(
+    object, predictions = NULL,
+    .levels = c(.5, .8, .95), ...,
+    .color_by = c("all_keys", "geo_value", "other_keys", ".response", "all", "none"),
+    .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
+    .base_color = "dodgerblue4",
+    .point_pred_color = "orange",
+    .max_facets = Inf) {
+  rlang::check_dots_empty()
+  arg_is_probabilities(.levels)
+  rlang::arg_match(.color_by)
+  rlang::arg_match(.facet_by)
+
+  if (!workflows::is_trained_workflow(object)) {
+    cli::cli_abort(c(
+      "Can't plot an untrained {.cls epi_workflow}.",
+      i = "Do you need to call `fit()`?"
+    ))
+  }
+
+  mold <- workflows::extract_mold(object)
+  y <- mold$outcomes
+  if (ncol(y) > 1) {
+    y <- y[, 1]
+    cli::cli_warn("Multiple outcome variables were detected. Displaying only 1.")
+  }
+  keys <- c("time_value", "geo_value", "key")
+  mold_roles <- names(mold$extras$roles)
+  edf <- dplyr::bind_cols(mold$extras$roles[mold_roles %in% keys], y)
+  if (starts_with_impl("ahead_", names(y))) {
+    old_name_y <- unlist(strsplit(names(y), "_"))
+    shift <- as.numeric(old_name_y[2])
+    new_name_y <- paste(old_name_y[-c(1:2)], collapse = "_")
+    edf <- dplyr::rename(edf, !!new_name_y := !!names(y))
+  } else if (starts_with_impl("lag_", names(y))) {
+    old_name_y <- unlist(strsplit(names(y), "_"))
+    shift <- -as.numeric(old_name_y[2])
+    new_name_y <- paste(old_name_y[-c(1:2)], collapse = "_")
+    edf <- dplyr::rename(edf, !!new_name_y := !!names(y))
+  }
+
+  if (!is.null(shift)) {
+    edf <- dplyr::mutate(edf, time_value = time_value + shift)
+  }
+  extra_keys <- setdiff(epi_keys_mold(mold), c("time_value", "geo_value"))
+  if (length(extra_keys) == 0L) extra_keys <- NULL
+  edf <- as_epi_df(edf,
+    as_of = object$fit$meta$as_of,
+    additional_metadata = list(other_keys = extra_keys)
+  )
+  if (is.null(predictions)) {
+    return(autoplot(
+      edf, new_name_y,
+      .color_by = .color_by, .facet_by = .facet_by, .base_color = .base_color,
+      .max_facets = .max_facets
+    ))
+  }
+
+  if ("target_date" %in% names(predictions)) {
+    if ("time_value" %in% names(predictions)) {
+      predictions <- dplyr::select(predictions, -time_value)
+    }
+    predictions <- dplyr::rename(predictions, time_value = target_date)
+  }
+  pred_cols_ok <- hardhat::check_column_names(predictions, epi_keys(edf))
+  if (!pred_cols_ok$ok) {
+    cli::cli_warn(c(
+      "`predictions` is missing required variables: {.var {pred_cols_ok$missing_names}}.",
+      i = "Plotting the original data."
+    ))
+    return(autoplot(
+      edf, !!new_name_y,
+      .color_by = .color_by, .facet_by = .facet_by, .base_color = .base_color,
+      .max_facets = .max_facets
+    ))
+  }
+
+  # First we plot the history, always faceted by everything
+  bp <- autoplot(edf, !!new_name_y,
+    .color_by = "none", .facet_by = "all_keys",
+    .base_color = "black", .max_facets = .max_facets
+  )
+
+  # Now, prepare matching facets in the predictions
+  ek <- kill_time_value(epi_keys(edf))
+  predictions <- predictions %>%
+    dplyr::mutate(
+      .facets = interaction(!!!rlang::syms(as.list(ek)), sep = "/"),
+    )
+  if (.max_facets < Inf) {
+    top_n <- levels(as.factor(bp$data$.facets))[seq_len(.max_facets)]
+    predictions <- dplyr::filter(predictions, .facets %in% top_n) %>%
+      dplyr::mutate(.facets = droplevels(.facets))
+  }
+
+
+  if (".pred_distn" %in% names(predictions)) {
+    bp <- plot_bands(bp, predictions, .levels, .base_color)
+  }
+
+  if (".pred" %in% names(predictions)) {
+    ntarget_dates <- dplyr::n_distinct(predictions$time_value)
+    if (ntarget_dates > 1L) {
+      bp <- bp +
+        ggplot2::geom_line(
+          data = predictions, ggplot2::aes(y = .data$.pred),
+          color = .point_pred_color
+        )
+    } else {
+      bp <- bp +
+        ggplot2::geom_point(
+          data = predictions, ggplot2::aes(y = .data$.pred),
+          color = .point_pred_color
+        )
+    }
+  }
+  bp
+}
+
+#' @export
+#' @rdname autoplot-epipred
+autoplot.canned_epipred <- function(
+    object, ...,
+    .color_by = c("all_keys", "geo_value", "other_keys", ".response", "all", "none"),
+    .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
+    .base_color = "dodgerblue4",
+    .point_pred_color = "orange",
+    .max_facets = Inf) {
+  rlang::check_dots_empty()
+  rlang::arg_match(.color_by)
+  rlang::arg_match(.facet_by)
+
+  ewf <- object$epi_workflow
+  predictions <- object$predictions %>%
+    dplyr::rename(time_value = target_date)
+
+  autoplot(ewf, predictions,
+    .color_by = .color_by, .facet_by = .facet_by,
+    .base_color = .base_color, .max_facets = .max_facets
+  )
+}
+
+starts_with_impl <- function(x, vars) {
+  n <- nchar(x)
+  x == substr(vars, 1, n)
+}
+
+plot_bands <- function(
+    base_plot, predictions,
+    levels = c(.5, .8, .95),
+    fill = "blue4",
+    alpha = 0.6,
+    linewidth = 0.05) {
+  innames <- names(predictions)
+  n <- length(levels)
+  alpha <- alpha / (n - 1)
+  l <- (1 - levels) / 2
+  l <- c(rev(l), 1 - l)
+
+  ntarget_dates <- dplyr::n_distinct(predictions$time_value)
+
+  predictions <- predictions %>%
+    dplyr::mutate(.pred_distn = dist_quantiles(quantile(.pred_distn, l), l)) %>%
+    pivot_quantiles_wider(.pred_distn)
+  qnames <- setdiff(names(predictions), innames)
+
+  for (i in 1:n) {
+    bottom <- qnames[i]
+    top <- rev(qnames)[i]
+    if (i == 1) {
+      if (ntarget_dates > 1L) {
+        base_plot <- base_plot +
+          ggplot2::geom_ribbon(
+            data = predictions,
+            ggplot2::aes(ymin = .data[[bottom]], ymax = .data[[top]]),
+            alpha = 0.2, linewidth = linewidth, fill = fill
+          )
+      } else {
+        base_plot <- base_plot +
+          ggplot2::geom_linerange(
+            data = predictions,
+            ggplot2::aes(ymin = .data[[bottom]], ymax = .data[[top]]),
+            alpha = 0.2, linewidth = 2, color = fill
+          )
+      }
+    } else {
+      if (ntarget_dates > 1L) {
+        base_plot <- base_plot +
+          ggplot2::geom_ribbon(
+            data = predictions,
+            ggplot2::aes(ymin = .data[[bottom]], ymax = .data[[top]]),
+            fill = fill, alpha = alpha
+          )
+      } else {
+        base_plot <- base_plot +
+          ggplot2::geom_linerange(
+            data = predictions,
+            ggplot2::aes(ymin = .data[[bottom]], ymax = .data[[top]]),
+            color = fill, alpha = alpha, linewidth = 2
+          )
+      }
+    }
+  }
+  base_plot
+}
+
+find_level <- function(x) {
+  unique((x < .5) * (1 - 2 * x) + (x > .5) * (1 - 2 * (1 - x)))
+}
diff --git a/R/epi_keys.R b/R/epi_keys.R
index 4a00cbd46..08e4595c3 100644
--- a/R/epi_keys.R
+++ b/R/epi_keys.R
@@ -24,7 +24,7 @@ epi_keys.data.frame <- function(x, other_keys = character(0L), ...) {
 
 #' @export
 epi_keys.epi_df <- function(x, ...) {
-  c("time_value", "geo_value", attributes(x)$metadata$other_keys)
+  c("time_value", "geo_value", attr(x, "metadata")$other_keys)
 }
 
 #' @export
@@ -32,6 +32,11 @@ epi_keys.recipe <- function(x, ...) {
   x$var_info$variable[x$var_info$role %in% c("time_value", "geo_value", "key")]
 }
 
+#' @export
+epi_keys.epi_workflow <- function(x, ...) {
+  epi_keys_mold(hardhat::extract_mold(x))
+}
+
 # a mold is a list extracted from a fitted workflow, gives info about
 # training data. But it doesn't have a class
 epi_keys_mold <- function(mold) {
@@ -45,3 +50,7 @@ kill_time_value <- function(v) {
   arg_is_chr(v)
   v[v != "time_value"]
 }
+
+epi_keys_only <- function(x, ...) {
+  kill_time_value(epi_keys(x, ...))
+}
diff --git a/R/workflow-printing.R b/R/workflow-printing.R
index c46d10848..d9c3446f9 100644
--- a/R/workflow-printing.R
+++ b/R/workflow-printing.R
@@ -105,7 +105,7 @@ print_preprocessor_formula <- function(x) {
   invisible(x)
 }
 
-print_prepocessor_variables <- function(x) {
+print_preprocessor_variables <- function(x) {
   variables <- workflows::extract_preprocessor(x)
   outcomes <- rlang::quo_get_expr(variables$outcomes)
   predictors <- rlang::quo_get_expr(variables$predictors)
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 07ceb6fec..c0537063a 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -47,9 +47,13 @@ reference:
   - title: Simple forecasters
     desc: Complete forecasters that produce reasonable baselines
     contents:
-    - contains("flatline")
-    - contains("arx")
-    - contains("cdc")
+    - contains("forecaster")
+    - contains("classifier")
+  - title: Forecaster modifications
+    desc: Constructors to modify forecaster arguments and utilities to produce `epi_workflow` objects
+    contents:
+    - contains("args_list")
+    - contains("_epi_workflow")
   - title: Helper functions for Hub submission
     contents:
     - flusight_hub_formatter
@@ -87,6 +91,10 @@ reference:
     contents:
     - contains("layer")
     - contains("slather")
+  - title: Automatic forecast visualization
+    contents:
+    - autoplot.epi_workflow
+    - autoplot.canned_epipred
   - title: Utilities for quantile distribution processing
     contents:
     - dist_quantiles
diff --git a/man/autoplot-epipred.Rd b/man/autoplot-epipred.Rd
new file mode 100644
index 000000000..80a552604
--- /dev/null
+++ b/man/autoplot-epipred.Rd
@@ -0,0 +1,123 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/autoplot.R
+\name{autoplot-epipred}
+\alias{autoplot-epipred}
+\alias{autoplot.epi_workflow}
+\alias{autoplot.canned_epipred}
+\title{Automatically plot an \code{epi_workflow} or \code{canned_epipred} object}
+\usage{
+\method{autoplot}{epi_workflow}(
+  object,
+  predictions = NULL,
+  .levels = c(0.5, 0.8, 0.95),
+  ...,
+  .color_by = c("all_keys", "geo_value", "other_keys", ".response", "all", "none"),
+  .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
+  .base_color = "dodgerblue4",
+  .point_pred_color = "orange",
+  .max_facets = Inf
+)
+
+\method{autoplot}{canned_epipred}(
+  object,
+  ...,
+  .color_by = c("all_keys", "geo_value", "other_keys", ".response", "all", "none"),
+  .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
+  .base_color = "dodgerblue4",
+  .point_pred_color = "orange",
+  .max_facets = Inf
+)
+}
+\arguments{
+\item{object}{An \code{epi_workflow}}
+
+\item{predictions}{A data frame with predictions. If \code{NULL}, only the
+original data is shown.}
+
+\item{.levels}{A numeric vector of levels to plot for any prediction bands.
+More than 3 levels begins to be difficult to see.}
+
+\item{...}{Ignored}
+
+\item{.color_by}{Which variables should determine the color(s) used to plot
+lines. Options include:
+\itemize{
+\item \code{all_keys} - the default uses the interaction of any key variables
+including the \code{geo_value}
+\item \code{geo_value} - \code{geo_value} only
+\item \code{other_keys} - any available keys that are not \code{geo_value}
+\item \code{.response} - the numeric variables (same as the y-axis)
+\item \code{all} - uses the interaction of all keys and numeric variables
+\item \code{none} - no coloring aesthetic is applied
+}}
+
+\item{.facet_by}{Similar to \code{.color_by} except that the default is to display
+each numeric variable on a separate facet}
+
+\item{.base_color}{If available, prediction bands will be shown with this
+color.}
+
+\item{.point_pred_color}{If available, point forecasts will be shown with this
+color.}
+
+\item{.max_facets}{Cut down of the number of facets displayed. Especially
+useful for testing when there are many \code{geo_value}'s or keys.}
+}
+\description{
+For a fit workflow, the training data will be displayed, the response by
+default. If \code{predictions} is not \code{NULL} then point and interval forecasts
+will be shown as well. Unfit workflows will result in an error, (you
+can simply call \code{autoplot()} on the original \code{epi_df}).
+}
+\examples{
+jhu <- case_death_rate_subset \%>\%
+  filter(time_value >= as.Date("2021-11-01"))
+
+r <- epi_recipe(jhu) \%>\%
+  step_epi_lag(death_rate, lag = c(0, 7, 14)) \%>\%
+  step_epi_ahead(death_rate, ahead = 7) \%>\%
+  step_epi_lag(case_rate, lag = c(0, 7, 14)) \%>\%
+  step_epi_naomit()
+
+f <- frosting() \%>\%
+  layer_residual_quantiles(
+    quantile_levels = c(.025, .1, .25, .75, .9, .975)
+  ) \%>\%
+  layer_threshold(dplyr::starts_with(".pred")) \%>\%
+  layer_add_target_date()
+
+wf <- epi_workflow(r, parsnip::linear_reg(), f) \%>\% fit(jhu)
+
+autoplot(wf)
+
+latest <- jhu \%>\% dplyr::filter(time_value >= max(time_value) - 14)
+preds <- predict(wf, latest)
+autoplot(wf, preds, .max_facets = 4)
+
+# ------- Show multiple horizons
+
+p <- lapply(c(7, 14, 21, 28), \(h) {
+  r <- epi_recipe(jhu) \%>\%
+    step_epi_lag(death_rate, lag = c(0, 7, 14)) \%>\%
+    step_epi_ahead(death_rate, ahead = h) \%>\%
+    step_epi_lag(case_rate, lag = c(0, 7, 14)) \%>\%
+    step_epi_naomit()
+  ewf <- epi_workflow(r, parsnip::linear_reg(), f) \%>\% fit(jhu)
+  td <- get_test_data(r, jhu)
+  predict(ewf, new_data = td)
+})
+
+p <- do.call(rbind, p)
+autoplot(wf, p, .max_facets = 4)
+
+# ------- Plotting canned forecaster output
+
+jhu <- case_death_rate_subset \%>\% filter(time_value >= as.Date("2021-11-01"))
+flat <- flatline_forecaster(jhu, "death_rate")
+autoplot(flat, .max_facets = 4)
+
+arx <- arx_forecaster(jhu, "death_rate", c("case_rate", "death_rate"),
+  args_list = arx_args_list(ahead = 14L)
+)
+autoplot(arx, .max_facets = 6)
+}
diff --git a/man/reexports.Rd b/man/reexports.Rd
index b23f00698..d376e1281 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -1,8 +1,9 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/reexports-tidymodels.R
+% Please edit documentation in R/autoplot.R, R/reexports-tidymodels.R
 \docType{import}
 \name{reexports}
 \alias{reexports}
+\alias{autoplot}
 \alias{fit}
 \alias{prep}
 \alias{bake}
@@ -15,6 +16,8 @@ below to see their documentation.
 \describe{
   \item{generics}{\code{\link[generics]{fit}}}
 
+  \item{ggplot2}{\code{\link[ggplot2]{autoplot}}}
+
   \item{recipes}{\code{\link[recipes]{bake}}, \code{\link[recipes]{prep}}}
 }}
 
diff --git a/tests/testthat/test-dist_quantiles.R b/tests/testthat/test-dist_quantiles.R
index 9c01ca103..ec885063b 100644
--- a/tests/testthat/test-dist_quantiles.R
+++ b/tests/testthat/test-dist_quantiles.R
@@ -26,6 +26,7 @@ test_that("single dist_quantiles works, quantiles are accessible", {
   )
 })
 
+
 test_that("quantile extrapolator works", {
   dstn <- dist_normal(c(10, 2), c(5, 10))
   qq <- extrapolate_quantiles(dstn, probs = c(.25, 0.5, .75))