diff --git a/DESCRIPTION b/DESCRIPTION
index 36a84f42..81a35b30 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: epipredict
 Title: Basic epidemiology forecasting methods
-Version: 0.1.14
+Version: 0.1.15
 Authors@R: c(
     person("Daniel J.", "McDonald", , "daniel@stat.ubc.ca", role = c("aut", "cre")),
     person("Ryan", "Tibshirani", , "ryantibs@cmu.edu", role = "aut"),
@@ -31,7 +31,7 @@ Imports:
     checkmate,
     cli,
     dplyr,
-    epiprocess (>= 0.10.4),
+    epiprocess (>= 0.11.2),
     generics,
     ggplot2,
     glue,
diff --git a/NAMESPACE b/NAMESPACE
index f39d9bfb..c2fa9494 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -47,6 +47,8 @@ S3method(forecast,epi_workflow)
 S3method(key_colnames,epi_workflow)
 S3method(key_colnames,recipe)
 S3method(mean,quantile_pred)
+S3method(plot,canned_epipred)
+S3method(plot,epi_workflow)
 S3method(predict,epi_workflow)
 S3method(predict,flatline)
 S3method(prep,check_enough_data)
diff --git a/NEWS.md b/NEWS.md
index 24c8e3b7..de698ee9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -23,6 +23,7 @@ Pre-1.0.0 numbering scheme: 0.x will indicate releases, while 0.0.x will indicat
 - Rename `check_enough_train_data()` to `check_enough_data()`, and generalize it
   enough to use as a check on either training or testing.
 - Add check for enough data to predict in `arx_forecaster()`
+- Adds the `.facet_filter` option in `epiprocess::autoplot()` (cmu-delphi/epiprocess#647).
 
 ## Improvements
 
diff --git a/R/autoplot.R b/R/autoplot.R
index 4e52ac3e..efad56ff 100644
--- a/R/autoplot.R
+++ b/R/autoplot.R
@@ -13,7 +13,7 @@ ggplot2::autoplot
 #'
 #'
 #' @inheritParams epiprocess::autoplot.epi_df
-#' @param object An `epi_workflow`
+#' @param object,x 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.
@@ -48,7 +48,7 @@ ggplot2::autoplot
 #'
 #' latest <- jhu %>% filter(time_value >= max(time_value) - 14)
 #' preds <- predict(wf, latest)
-#' autoplot(wf, preds, .max_facets = 4)
+#' autoplot(wf, preds, .facet_filter = geo_value %in% c("ca", "ny", "de", "mt"))
 #'
 #' # ------- Show multiple horizons
 #'
@@ -63,19 +63,19 @@ ggplot2::autoplot
 #' })
 #'
 #' p <- do.call(rbind, p)
-#' autoplot(wf, p, .max_facets = 4)
+#' autoplot(wf, p, .facet_filter = geo_value %in% c("ca", "ny", "de", "mt"))
 #'
 #' # ------- Plotting canned forecaster output
 #'
 #' jhu <- covid_case_death_rates %>%
 #'   filter(time_value >= as.Date("2021-11-01"))
 #' flat <- flatline_forecaster(jhu, "death_rate")
-#' autoplot(flat, .max_facets = 4)
+#' autoplot(flat, .facet_filter = geo_value %in% c("ca", "ny", "de", "mt"))
 #'
 #' arx <- arx_forecaster(jhu, "death_rate", c("case_rate", "death_rate"),
 #'   args_list = arx_args_list(ahead = 14L)
 #' )
-#' autoplot(arx, .max_facets = 6)
+#' autoplot(arx, .facet_filter = geo_value %in% c("ca", "ny", "de", "mt", "mo", "in"))
 NULL
 
 #' @export
@@ -87,7 +87,8 @@ autoplot.epi_workflow <- function(
     .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
     .base_color = "dodgerblue4",
     .point_pred_color = "orange",
-    .max_facets = Inf) {
+    .facet_filter = NULL,
+    .max_facets = deprecated()) {
   rlang::check_dots_empty()
   arg_is_probabilities(.levels)
   rlang::arg_match(.color_by)
@@ -134,7 +135,7 @@ autoplot.epi_workflow <- function(
     return(autoplot(
       edf, new_name_y,
       .color_by = .color_by, .facet_by = .facet_by, .base_color = .base_color,
-      .max_facets = .max_facets
+      .facet_filter = {{ .facet_filter }}
     ))
   }
 
@@ -153,25 +154,25 @@ autoplot.epi_workflow <- function(
     return(autoplot(
       edf, !!new_name_y,
       .color_by = .color_by, .facet_by = .facet_by, .base_color = .base_color,
-      .max_facets = .max_facets
+      .facet_filter = {{ .facet_filter }}
     ))
   }
 
   # 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
+    .base_color = "black", .facet_filter =  {{ .facet_filter }}
   )
 
   # Now, prepare matching facets in the predictions
   ek <- epi_keys_only(edf)
   predictions <- predictions %>%
     mutate(
-      .facets = interaction(!!!rlang::syms(as.list(ek)), sep = "/"),
+      .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 <- filter(predictions, .facets %in% top_n) %>%
+  .facet_filter <- rlang::enquo(.facet_filter)
+  if (!rlang::quo_is_null(.facet_filter) && ".facets" %in% names(bp$data)) {
+    predictions <- filter(predictions, .facets %in% unique(bp$data$.facets)) %>%
       mutate(.facets = droplevels(.facets))
   }
 
@@ -207,7 +208,8 @@ autoplot.canned_epipred <- function(
     .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
     .base_color = "dodgerblue4",
     .point_pred_color = "orange",
-    .max_facets = Inf) {
+    .facet_filter = NULL,
+    .max_facets = deprecated()) {
   rlang::check_dots_empty()
   rlang::arg_match(.color_by)
   rlang::arg_match(.facet_by)
@@ -218,10 +220,23 @@ autoplot.canned_epipred <- function(
 
   autoplot(ewf, predictions,
     .color_by = .color_by, .facet_by = .facet_by,
-    .base_color = .base_color, .max_facets = .max_facets
+    .base_color = .base_color, .facet_filter = {{ .facet_filter }}
   )
 }
 
+#' @export
+#' @rdname autoplot-epipred
+plot.epi_workflow <- function(x, ...) {
+  autoplot(x, ...)
+}
+
+#' @export
+#' @rdname autoplot-epipred
+plot.canned_epipred <- function(x, ...) {
+  autoplot(x, ...)
+}
+
+
 starts_with_impl <- function(x, vars) {
   n <- nchar(x)
   x == substr(vars, 1, n)
diff --git a/R/step_climate.R b/R/step_climate.R
index 2ae2c79d..6e2817fa 100644
--- a/R/step_climate.R
+++ b/R/step_climate.R
@@ -266,7 +266,7 @@ prep.step_climate <- function(x, training, info = NULL, ...) {
 
   fn <- switch(x$center_method,
     mean = function(x, w) stats::weighted.mean(x, w, na.rm = TRUE),
-    median = function(x, w) median(x, na.rm = TRUE)
+    median = function(x, w) stats::median(x, na.rm = TRUE)
   )
   # suppose it's week 52, and there is no week 53 this year; then
   # as originally written for 1 week ahead this grabs from week 52+1 %% 53
diff --git a/man/autoplot-epipred.Rd b/man/autoplot-epipred.Rd
index 087f25db..f4657967 100644
--- a/man/autoplot-epipred.Rd
+++ b/man/autoplot-epipred.Rd
@@ -4,6 +4,8 @@
 \alias{autoplot-epipred}
 \alias{autoplot.epi_workflow}
 \alias{autoplot.canned_epipred}
+\alias{plot.epi_workflow}
+\alias{plot.canned_epipred}
 \title{Automatically plot an \code{epi_workflow} or \code{canned_epipred} object}
 \usage{
 \method{autoplot}{epi_workflow}(
@@ -15,7 +17,8 @@
   .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
   .base_color = "dodgerblue4",
   .point_pred_color = "orange",
-  .max_facets = Inf
+  .facet_filter = NULL,
+  .max_facets = deprecated()
 )
 
 \method{autoplot}{canned_epipred}(
@@ -25,11 +28,16 @@
   .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
   .base_color = "dodgerblue4",
   .point_pred_color = "orange",
-  .max_facets = Inf
+  .facet_filter = NULL,
+  .max_facets = deprecated()
 )
+
+\method{plot}{epi_workflow}(x, ...)
+
+\method{plot}{canned_epipred}(x, ...)
 }
 \arguments{
-\item{object}{An \code{epi_workflow}}
+\item{object, x}{An \code{epi_workflow}}
 
 \item{predictions}{A data frame with predictions. If \code{NULL}, only the
 original data is shown.}
@@ -60,8 +68,16 @@ 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.}
+\item{.facet_filter}{Select which facets will be displayed. Especially
+useful for when there are many \code{geo_value}'s or keys. This is a
+<\code{\link[=args_data_masking]{rlang}}> expression along the lines of \code{\link[dplyr:filter]{dplyr::filter()}}.
+However, it must be a single expression combined with the \code{&} operator. This
+contrasts to the typical use case which allows multiple comma-separated expressions
+which are implicitly combined with \code{&}. When multiple variables are selected
+with \code{...}, their names can be filtered in combination with other factors
+by using \code{.response_name}. See the examples below.}
+
+\item{.max_facets}{\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}}}
 }
 \description{
 For a fit workflow, the training data will be displayed, the response by
@@ -90,7 +106,7 @@ autoplot(wf)
 
 latest <- jhu \%>\% filter(time_value >= max(time_value) - 14)
 preds <- predict(wf, latest)
-autoplot(wf, preds, .max_facets = 4)
+autoplot(wf, preds, .facet_filter = geo_value \%in\% c("ca", "ny", "de", "mt"))
 
 # ------- Show multiple horizons
 
@@ -105,17 +121,17 @@ p <- lapply(c(7, 14, 21, 28), function(h) {
 })
 
 p <- do.call(rbind, p)
-autoplot(wf, p, .max_facets = 4)
+autoplot(wf, p, .facet_filter = geo_value \%in\% c("ca", "ny", "de", "mt"))
 
 # ------- Plotting canned forecaster output
 
 jhu <- covid_case_death_rates \%>\%
   filter(time_value >= as.Date("2021-11-01"))
 flat <- flatline_forecaster(jhu, "death_rate")
-autoplot(flat, .max_facets = 4)
+autoplot(flat, .facet_filter = geo_value \%in\% c("ca", "ny", "de", "mt"))
 
 arx <- arx_forecaster(jhu, "death_rate", c("case_rate", "death_rate"),
   args_list = arx_args_list(ahead = 14L)
 )
-autoplot(arx, .max_facets = 6)
+autoplot(arx, .facet_filter = geo_value \%in\% c("ca", "ny", "de", "mt", "mo", "in"))
 }