diff --git a/NAMESPACE b/NAMESPACE
index d5a0a8a9de..efeddbb12d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -165,11 +165,13 @@ export(GeomBlank)
 export(GeomBoxplot)
 export(GeomCol)
 export(GeomContour)
+export(GeomContourFilled)
 export(GeomCrossbar)
 export(GeomCurve)
 export(GeomCustomAnn)
 export(GeomDensity)
 export(GeomDensity2d)
+export(GeomDensity2dFilled)
 export(GeomDotplot)
 export(GeomErrorbar)
 export(GeomErrorbarh)
@@ -231,6 +233,7 @@ export(StatContourFilled)
 export(StatCount)
 export(StatDensity)
 export(StatDensity2d)
+export(StatDensity2dFilled)
 export(StatEcdf)
 export(StatEllipse)
 export(StatFunction)
@@ -341,7 +344,9 @@ export(geom_crossbar)
 export(geom_curve)
 export(geom_density)
 export(geom_density2d)
+export(geom_density2d_filled)
 export(geom_density_2d)
+export(geom_density_2d_filled)
 export(geom_dotplot)
 export(geom_errorbar)
 export(geom_errorbarh)
@@ -591,7 +596,9 @@ export(stat_contour_filled)
 export(stat_count)
 export(stat_density)
 export(stat_density2d)
+export(stat_density2d_filled)
 export(stat_density_2d)
+export(stat_density_2d_filled)
 export(stat_ecdf)
 export(stat_ellipse)
 export(stat_function)
diff --git a/NEWS.md b/NEWS.md
index 8e29f6bf5e..32578aa3ad 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,6 +3,10 @@
 * `annotation_raster()` adds support for native rasters. For large rasters,
   native rasters render significantly faster than arrays (@kent37, #3388)
   
+* A newly added geom `geom_density_2d_filled()` and associated stat 
+  `stat_density_2d_filled()` can draw filled density contours
+  (@clauswilke, #3846).
+  
 * Support graphics devices that use the `file` argument instead of `fileneame` 
   in `ggsave()` (@bwiernik, #3810)
 
diff --git a/R/geom-contour.r b/R/geom-contour.r
index 2dc81cba4a..24a79c938b 100644
--- a/R/geom-contour.r
+++ b/R/geom-contour.r
@@ -1,14 +1,15 @@
-#' 2d contours of a 3d surface
+#' 2D contours of a 3D surface
 #'
-#' ggplot2 can not draw true 3d surfaces, but you can use `geom_contour`
-#' and [geom_tile()] to visualise 3d surfaces in 2d. To be a valid
-#' surface, the data must contain only a single row for each unique combination
-#' of the variables mapped to the `x` and `y` aesthetics. Contouring
+#' ggplot2 can not draw true 3D surfaces, but you can use `geom_contour()`,
+#' `geom_contour_filled()`, and [geom_tile()] to visualise 3D surfaces in 2D.
+#' To specify a valid surface, the data must contain `x`, `y`, and `z` coordinates,
+#' and each unique combination of `x` and `y` can appear exactly once. Contouring
 #' tends to work best when `x` and `y` form a (roughly) evenly
 #' spaced grid. If your data is not evenly spaced, you may want to interpolate
-#' to a grid before visualising.
+#' to a grid before visualising, see [geom_density_2d()].
 #'
 #' @eval rd_aesthetics("geom", "contour")
+#' @eval rd_aesthetics("geom", "contour_filled")
 #' @inheritParams layer
 #' @inheritParams geom_point
 #' @inheritParams geom_path
@@ -20,7 +21,7 @@
 #' @seealso [geom_density_2d()]: 2d density contours
 #' @export
 #' @examples
-#' #' # Basic plot
+#' # Basic plot
 #' v <- ggplot(faithfuld, aes(waiting, eruptions, z = density))
 #' v + geom_contour()
 #'
@@ -33,7 +34,7 @@
 #' v + geom_contour_filled()
 #'
 #' # Setting bins creates evenly spaced contours in the range of the data
-#' v + geom_contour(bins = 2)
+#' v + geom_contour(bins = 5)
 #' v + geom_contour(bins = 10)
 #'
 #' # Setting binwidth does the same thing, parameterised by the distance
@@ -95,7 +96,7 @@ geom_contour_filled <- function(mapping = NULL, data = NULL,
     data = data,
     mapping = mapping,
     stat = stat,
-    geom = GeomPolygon,
+    geom = GeomContourFilled,
     position = position,
     show.legend = show.legend,
     inherit.aes = inherit.aes,
@@ -123,3 +124,11 @@ GeomContour <- ggproto("GeomContour", GeomPath,
     alpha = NA
   )
 )
+
+#' @rdname ggplot2-ggproto
+#' @format NULL
+#' @usage NULL
+#' @export
+#' @include geom-polygon.r
+GeomContourFilled <- ggproto("GeomContourFilled", GeomPolygon)
+
diff --git a/R/geom-density2d.r b/R/geom-density2d.r
index 2fb8cb5967..b7dbb7a474 100644
--- a/R/geom-density2d.r
+++ b/R/geom-density2d.r
@@ -1,27 +1,41 @@
-#' Contours of a 2d density estimate
+#' Contours of a 2D density estimate
 #'
 #' Perform a 2D kernel density estimation using [MASS::kde2d()] and
 #' display the results with contours. This can be useful for dealing with
-#' overplotting. This is a 2d version of [geom_density()].
+#' overplotting. This is a 2D version of [geom_density()]. `geom_density_2d()`
+#' draws contour lines, and `geom_density_2d_filled()` draws filled contour
+#' bands.
 #'
 #' @eval rd_aesthetics("geom", "density_2d")
-#' @seealso [geom_contour()] for information about how contours
-#'  are drawn; [geom_bin2d()] for another way of dealing with
+#' @eval rd_aesthetics("geom", "density_2d_filled")
+#' @seealso [geom_contour()], [geom_contour_filled()] for information about
+#'  how contours are drawn; [geom_bin2d()] for another way of dealing with
 #'  overplotting.
 #' @param geom,stat Use to override the default connection between
 #'   `geom_density_2d` and `stat_density_2d`.
 #' @inheritParams layer
 #' @inheritParams geom_point
 #' @inheritParams geom_path
+#' @param contour_var Character string identifying the variable to contour
+#'   by. Can be one of `"density"`, `"ndensity"`, or `"count"`. See the section
+#'   on computed variables for details.
 #' @export
 #' @examples
 #' m <- ggplot(faithful, aes(x = eruptions, y = waiting)) +
 #'  geom_point() +
 #'  xlim(0.5, 6) +
 #'  ylim(40, 110)
+#'
+#' # contour lines
 #' m + geom_density_2d()
+#'
 #' \donttest{
-#' m + stat_density_2d(aes(fill = after_stat(level)), geom = "polygon")
+#' # contour bands
+#' m + geom_density_2d_filled(alpha = 0.5)
+#'
+#' # contour bands and contour lines
+#' m + geom_density_2d_filled(alpha = 0.5) +
+#'   geom_density_2d(size = 0.25, colour = "black")
 #'
 #' set.seed(4393)
 #' dsmall <- diamonds[sample(nrow(diamonds), 1000), ]
@@ -30,23 +44,29 @@
 #' # set of contours for each value of that variable
 #' d + geom_density_2d(aes(colour = cut))
 #'
-#' # Similarly, if you apply faceting to the plot, contours will be
-#' # drawn for each facet, but the levels will calculated across all facets
-#' d + stat_density_2d(aes(fill = after_stat(level)), geom = "polygon") +
-#'   facet_grid(. ~ cut) + scale_fill_viridis_c()
-#' # To override this behavior (for instace, to better visualize the density
-#' # within each facet), use after_stat(nlevel)
-#' d + stat_density_2d(aes(fill = after_stat(nlevel)), geom = "polygon") +
-#'   facet_grid(. ~ cut) + scale_fill_viridis_c()
+#' # If you draw filled contours across multiple facets, the same bins are
+#' # used across all facets
+#' d + geom_density_2d_filled() + facet_wrap(vars(cut))
+#' # If you want to make sure the peak intensity is the same in each facet,
+#' # use `contour_var = "ndensity"`.
+#' d + geom_density_2d_filled(contour_var = "ndensity") + facet_wrap(vars(cut))
+#' # If you want to scale intensity by the number of observations in each group,
+#' # use `contour_var = "count"`.
+#' d + geom_density_2d_filled(contour_var = "count") + facet_wrap(vars(cut))
 #'
-#' # If we turn contouring off, we can use use geoms like tiles:
-#' d + stat_density_2d(geom = "raster", aes(fill = after_stat(density)), contour = FALSE)
+#' # If we turn contouring off, we can use other geoms, such as tiles:
+#' d + stat_density_2d(
+#'   geom = "raster",
+#'   aes(fill = after_stat(density)),
+#'   contour = FALSE
+#' ) + scale_fill_viridis_c()
 #' # Or points:
 #' d + stat_density_2d(geom = "point", aes(size = after_stat(density)), n = 20, contour = FALSE)
 #' }
 geom_density_2d <- function(mapping = NULL, data = NULL,
-                            stat = "density2d", position = "identity",
+                            stat = "density_2d", position = "identity",
                             ...,
+                            contour_var = "density",
                             lineend = "butt",
                             linejoin = "round",
                             linemitre = 10,
@@ -65,6 +85,8 @@ geom_density_2d <- function(mapping = NULL, data = NULL,
       lineend = lineend,
       linejoin = linejoin,
       linemitre = linemitre,
+      contour = TRUE,
+      contour_var = contour_var,
       na.rm = na.rm,
       ...
     )
@@ -84,3 +106,42 @@ geom_density2d <- geom_density_2d
 GeomDensity2d <- ggproto("GeomDensity2d", GeomPath,
   default_aes = aes(colour = "#3366FF", size = 0.5, linetype = 1, alpha = NA)
 )
+
+#' @export
+#' @rdname geom_density_2d
+geom_density_2d_filled <- function(mapping = NULL, data = NULL,
+                                   stat = "density_2d_filled", position = "identity",
+                                   ...,
+                                   contour_var = "density",
+                                   na.rm = FALSE,
+                                   show.legend = NA,
+                                   inherit.aes = TRUE) {
+  layer(
+    data = data,
+    mapping = mapping,
+    stat = stat,
+    geom = GeomDensity2dFilled,
+    position = position,
+    show.legend = show.legend,
+    inherit.aes = inherit.aes,
+    params = list(
+      na.rm = na.rm,
+      contour = TRUE,
+      contour_var = contour_var,
+      ...
+    )
+  )
+}
+
+#' @export
+#' @rdname geom_density_2d
+#' @usage NULL
+geom_density2d_filled <- geom_density_2d_filled
+
+#' @rdname ggplot2-ggproto
+#' @format NULL
+#' @usage NULL
+#' @export
+#' @include geom-polygon.r
+GeomDensity2dFilled <- ggproto("GeomDensity2dFilled", GeomPolygon)
+
diff --git a/R/stat-contour.r b/R/stat-contour.r
index 25f96bfe88..5ef4dedc4b 100644
--- a/R/stat-contour.r
+++ b/R/stat-contour.r
@@ -2,11 +2,21 @@
 #' @inheritParams geom_contour
 #' @export
 #' @eval rd_aesthetics("stat", "contour")
+#' @eval rd_aesthetics("stat", "contour_filled")
 #' @section Computed variables:
+#' The computed variables differ somewhat for contour lines (computed by
+#' `stat_contour()`) and contour bands (filled contours, computed by `stat_contour_filled()`).
+#' The variables `nlevel` and `piece` are available for both, whereas `level_low`, `level_high`,
+#' and `level_mid` are only available for bands. The variable `level` is a numeric or a factor
+#' depending on whether lines or bands are calculated.
 #' \describe{
-#'  \item{level}{height of contour}
-#'  \item{nlevel}{height of contour, scaled to maximum of 1}
-#'  \item{piece}{contour piece (an integer)}
+#'  \item{`level`}{Height of contour. For contour lines, this is numeric vector that
+#'    represents bin boundaries. For contour bands, this is an ordered factor that
+#'    represents bin ranges.}
+#'  \item{`level_low`, `level_high`, `level_mid`}{(contour bands only) Lower and upper
+#'    bin boundaries for each band, as well the mid point between the boundaries.}
+#'  \item{`nlevel`}{Height of contour, scaled to maximum of 1.}
+#'  \item{`piece`}{Contour piece (an integer).}
 #' }
 #' @rdname geom_contour
 stat_contour <- function(mapping = NULL, data = NULL,
@@ -39,7 +49,7 @@ stat_contour <- function(mapping = NULL, data = NULL,
 #' @rdname geom_contour
 #' @export
 stat_contour_filled <- function(mapping = NULL, data = NULL,
-                                geom = "polygon", position = "identity",
+                                geom = "contour_filled", position = "identity",
                                 ...,
                                 bins = NULL,
                                 binwidth = NULL,
@@ -74,11 +84,15 @@ StatContour <- ggproto("StatContour", Stat,
   required_aes = c("x", "y", "z"),
   default_aes = aes(order = after_stat(level)),
 
-  compute_group = function(data, scales, bins = NULL, binwidth = NULL,
+  setup_params = function(data, params) {
+    params$z.range <- range(data$z, na.rm = TRUE, finite = TRUE)
+    params
+  },
+
+  compute_group = function(data, scales, z.range, bins = NULL, binwidth = NULL,
                            breaks = NULL, na.rm = FALSE) {
 
-    z_range <- range(data$z, na.rm = TRUE, finite = TRUE)
-    breaks <- contour_breaks(z_range, bins, binwidth, breaks)
+    breaks <- contour_breaks(z.range, bins, binwidth, breaks)
 
     isolines <- xyz_to_isolines(data, breaks)
     path_df <- iso_to_path(isolines, data$group[1])
@@ -99,16 +113,23 @@ StatContourFilled <- ggproto("StatContourFilled", Stat,
   required_aes = c("x", "y", "z"),
   default_aes = aes(order = after_stat(level), fill = after_stat(level)),
 
-  compute_group = function(data, scales, bins = NULL, binwidth = NULL, breaks = NULL, na.rm = FALSE) {
+  setup_params = function(data, params) {
+    params$z.range <- range(data$z, na.rm = TRUE, finite = TRUE)
+    params
+  },
 
-    z_range <- range(data$z, na.rm = TRUE, finite = TRUE)
-    breaks <- contour_breaks(z_range, bins, binwidth, breaks)
+  compute_group = function(data, scales, z.range, bins = NULL, binwidth = NULL, breaks = NULL, na.rm = FALSE) {
+    breaks <- contour_breaks(z.range, bins, binwidth, breaks)
 
     isobands <- xyz_to_isobands(data, breaks)
     names(isobands) <- pretty_isoband_levels(names(isobands))
     path_df <- iso_to_polygon(isobands, data$group[1])
 
     path_df$level <- ordered(path_df$level, levels = names(isobands))
+    path_df$level_low <- breaks[as.numeric(path_df$level)]
+    path_df$level_high <- breaks[as.numeric(path_df$level) + 1]
+    path_df$level_mid <- 0.5*(path_df$level_low + path_df$level_high)
+    path_df$nlevel <- rescale_max(path_df$level_high)
 
     path_df
   }
diff --git a/R/stat-density-2d.r b/R/stat-density-2d.r
index 1a5fcf0e7a..85c8118376 100644
--- a/R/stat-density-2d.r
+++ b/R/stat-density-2d.r
@@ -1,8 +1,11 @@
 #' @export
 #' @rdname geom_density_2d
 #' @param contour If `TRUE`, contour the results of the 2d density
-#'   estimation
-#' @param n number of grid points in each direction
+#'   estimation.
+#' @param contour_var Character string identifying the variable to contour
+#'   by. Can be one of `"density"`, `"ndensity"`, or `"count"`. See the section
+#'   on computed variables for details.
+#' @param n Number of grid points in each direction.
 #' @param h Bandwidth (vector of length two). If `NULL`, estimated
 #'   using [MASS::bandwidth.nrd()].
 #' @param adjust A multiplicative bandwidth adjustment to be used if 'h' is
@@ -10,17 +13,29 @@
 #'    using the a bandwidth estimator. For example, `adjust = 1/2` means
 #'    use half of the default bandwidth.
 #' @section Computed variables:
-#' Same as [stat_contour()]
-#'
-#' With the addition of:
+#' `stat_density_2d()` and `stat_density_2d_filled()` compute different
+#' variables depending on whether contouring is turned on or off. With
+#' contouring off (`contour = FALSE`), both stats behave the same, and the
+#' following variables are provided:
 #' \describe{
-#'   \item{density}{the density estimate}
-#'   \item{ndensity}{density estimate, scaled to maximum of 1}
+#'   \item{`density`}{The density estimate.}
+#'   \item{`ndensity`}{Density estimate, scaled to a maximum of 1.}
+#'   \item{`count`}{Density estimate * number of observations in group.}
+#'   \item{`n`}{Number of observations in each group.}
 #' }
+#'
+#' With contouring on (`contour = TRUE`), either [stat_contour()] or
+#' [stat_contour_filled()] (for contour lines or contour bands,
+#' respectively) is run after the density estimate has been obtained,
+#' and the computed variables are determined by these stats.
+#' Contours are calculated for one of the three types of density estimates
+#' obtained before contouring, `density`, `ndensity`, and `count`. Which
+#' of those should be used is determined by the `contour_var` parameter.
 stat_density_2d <- function(mapping = NULL, data = NULL,
                             geom = "density_2d", position = "identity",
                             ...,
                             contour = TRUE,
+                            contour_var = "density",
                             n = 100,
                             h = NULL,
                             adjust = c(1, 1),
@@ -38,6 +53,7 @@ stat_density_2d <- function(mapping = NULL, data = NULL,
     params = list(
       na.rm = na.rm,
       contour = contour,
+      contour_var = contour_var,
       n = n,
       h = h,
       adjust = adjust,
@@ -46,11 +62,50 @@ stat_density_2d <- function(mapping = NULL, data = NULL,
   )
 }
 
-#' @export
 #' @rdname geom_density_2d
 #' @usage NULL
+#' @export
 stat_density2d <- stat_density_2d
 
+#' @rdname geom_density_2d
+#' @export
+stat_density_2d_filled <- function(mapping = NULL, data = NULL,
+                                   geom = "density_2d_filled", position = "identity",
+                                   ...,
+                                   contour = TRUE,
+                                   contour_var = "density",
+                                   n = 100,
+                                   h = NULL,
+                                   adjust = c(1, 1),
+                                   na.rm = FALSE,
+                                   show.legend = NA,
+                                   inherit.aes = TRUE) {
+  layer(
+    data = data,
+    mapping = mapping,
+    stat = StatDensity2dFilled,
+    geom = geom,
+    position = position,
+    show.legend = show.legend,
+    inherit.aes = inherit.aes,
+    params = list(
+      na.rm = na.rm,
+      contour = contour,
+      contour_var = contour_var,
+      n = n,
+      h = h,
+      adjust = adjust,
+      ...
+    )
+  )
+}
+
+#' @rdname geom_density_2d
+#' @usage NULL
+#' @export
+stat_density2d_filled <- stat_density_2d_filled
+
+
 #' @rdname ggplot2-ggproto
 #' @format NULL
 #' @usage NULL
@@ -60,30 +115,85 @@ StatDensity2d <- ggproto("StatDensity2d", Stat,
 
   required_aes = c("x", "y"),
 
+  extra_params = c(
+    "na.rm", "contour", "contour_var",
+    "bins", "binwidth", "breaks"
+  ),
+
+  # when contouring is on, are we returning lines or bands?
+  contour_type = "lines",
+
+  compute_layer = function(self, data, params, layout) {
+    # first run the regular layer calculation to infer densities
+    data <- ggproto_parent(Stat, self)$compute_layer(data, params, layout)
+
+    # if we're not contouring we're done
+    if (!isTRUE(params$contour)) return(data)
+
+    # set up data and parameters for contouring
+    contour_var <- params$contour_var %||% "density"
+    if (!isTRUE(contour_var %in% c("density", "ndensity", "count"))) {
+      abort(glue(
+        'Unsupported value for `contour_var`: {contour_var}\n',
+        'Supported values are "density", "ndensity", and "count".'
+      ))
+    }
+    data$z <- data[[contour_var]]
+    z.range <- range(data$z, na.rm = TRUE, finite = TRUE)
+    params <- params[intersect(names(params), c("bins", "binwidth", "breaks"))]
+    params$z.range <- z.range
+
+    if (isTRUE(self$contour_type == "bands")) {
+      contour_stat <- StatContourFilled
+    } else { # lines is the default
+      contour_stat <- StatContour
+    }
+
+    args <- c(list(data = quote(data), scales = quote(scales)), params)
+    dapply(data, "PANEL", function(data) {
+      scales <- layout$get_scales(data$PANEL[1])
+      tryCatch(do.call(contour_stat$compute_panel, args), error = function(e) {
+        warn(glue("Computation failed in `{snake_class(self)}()`:\n{e$message}"))
+        new_data_frame()
+      })
+    })
+  },
+
   compute_group = function(data, scales, na.rm = FALSE, h = NULL, adjust = c(1, 1),
-                           contour = TRUE, n = 100, bins = NULL,
-                           binwidth = NULL) {
+                           n = 100, ...) {
     if (is.null(h)) {
       h <- c(MASS::bandwidth.nrd(data$x), MASS::bandwidth.nrd(data$y))
       h <- h * adjust
     }
 
+    # calculate density
     dens <- MASS::kde2d(
       data$x, data$y, h = h, n = n,
       lims = c(scales$x$dimension(), scales$y$dimension())
     )
+
+    # prepare final output data frame
+    nx <- nrow(data) # number of observations in this group
     df <- expand.grid(x = dens$x, y = dens$y)
-    df$z <- as.vector(dens$z)
+    df$density <- as.vector(dens$z)
     df$group <- data$group[1]
-
-    if (contour) {
-        StatContour$compute_panel(df, scales, bins, binwidth)
-    } else {
-      names(df) <- c("x", "y", "density", "group")
-      df$ndensity <- df$density / max(df$density, na.rm = TRUE)
-      df$level <- 1
-      df$piece <- 1
-      df
-    }
+    df$ndensity <- df$density / max(df$density, na.rm = TRUE)
+    df$count <- nx * df$density
+    df$n <- nx
+    df$level <- 1
+    df$piece <- 1
+    df
   }
 )
+
+
+
+#' @rdname ggplot2-ggproto
+#' @format NULL
+#' @usage NULL
+#' @export
+StatDensity2dFilled <- ggproto("StatDensity2dFilled", StatDensity2d,
+  default_aes = aes(colour = NA, fill = after_stat(level)),
+  contour_type = "bands"
+)
+
diff --git a/man/geom_contour.Rd b/man/geom_contour.Rd
index 34c88bdd64..e79df6c41d 100644
--- a/man/geom_contour.Rd
+++ b/man/geom_contour.Rd
@@ -5,7 +5,7 @@
 \alias{geom_contour_filled}
 \alias{stat_contour}
 \alias{stat_contour_filled}
-\title{2d contours of a 3d surface}
+\title{2D contours of a 3D surface}
 \usage{
 geom_contour(
   mapping = NULL,
@@ -55,7 +55,7 @@ stat_contour(
 stat_contour_filled(
   mapping = NULL,
   data = NULL,
-  geom = "polygon",
+  geom = "contour_filled",
   position = "identity",
   ...,
   bins = NULL,
@@ -129,13 +129,13 @@ the default plot specification, e.g. \code{\link[=borders]{borders()}}.}
 \item{geom}{The geometric object to use display the data}
 }
 \description{
-ggplot2 can not draw true 3d surfaces, but you can use \code{geom_contour}
-and \code{\link[=geom_tile]{geom_tile()}} to visualise 3d surfaces in 2d. To be a valid
-surface, the data must contain only a single row for each unique combination
-of the variables mapped to the \code{x} and \code{y} aesthetics. Contouring
+ggplot2 can not draw true 3D surfaces, but you can use \code{geom_contour()},
+\code{geom_contour_filled()}, and \code{\link[=geom_tile]{geom_tile()}} to visualise 3D surfaces in 2D.
+To specify a valid surface, the data must contain \code{x}, \code{y}, and \code{z} coordinates,
+and each unique combination of \code{x} and \code{y} can appear exactly once. Contouring
 tends to work best when \code{x} and \code{y} form a (roughly) evenly
 spaced grid. If your data is not evenly spaced, you may want to interpolate
-to a grid before visualising.
+to a grid before visualising, see \code{\link[=geom_density_2d]{geom_density_2d()}}.
 }
 \section{Aesthetics}{
 
@@ -153,6 +153,21 @@ to a grid before visualising.
 Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
 
 
+\code{geom_contour_filled()} understands the following aesthetics (required aesthetics are in bold):
+\itemize{
+\item \strong{\code{x}}
+\item \strong{\code{y}}
+\item \code{alpha}
+\item \code{colour}
+\item \code{fill}
+\item \code{group}
+\item \code{linetype}
+\item \code{size}
+\item \code{subgroup}
+}
+Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
+
+
 \code{stat_contour()} understands the following aesthetics (required aesthetics are in bold):
 \itemize{
 \item \strong{\code{x}}
@@ -162,19 +177,40 @@ Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
 \item \code{order}
 }
 Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
+
+
+\code{stat_contour_filled()} understands the following aesthetics (required aesthetics are in bold):
+\itemize{
+\item \strong{\code{x}}
+\item \strong{\code{y}}
+\item \strong{\code{z}}
+\item \code{fill}
+\item \code{group}
+\item \code{order}
+}
+Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
 }
 
 \section{Computed variables}{
 
+The computed variables differ somewhat for contour lines (computed by
+\code{stat_contour()}) and contour bands (filled contours, computed by \code{stat_contour_filled()}).
+The variables \code{nlevel} and \code{piece} are available for both, whereas \code{level_low}, \code{level_high},
+and \code{level_mid} are only available for bands. The variable \code{level} is a numeric or a factor
+depending on whether lines or bands are calculated.
 \describe{
-\item{level}{height of contour}
-\item{nlevel}{height of contour, scaled to maximum of 1}
-\item{piece}{contour piece (an integer)}
+\item{\code{level}}{Height of contour. For contour lines, this is numeric vector that
+represents bin boundaries. For contour bands, this is an ordered factor that
+represents bin ranges.}
+\item{\code{level_low}, \code{level_high}, \code{level_mid}}{(contour bands only) Lower and upper
+bin boundaries for each band, as well the mid point between the boundaries.}
+\item{\code{nlevel}}{Height of contour, scaled to maximum of 1.}
+\item{\code{piece}}{Contour piece (an integer).}
 }
 }
 
 \examples{
-#' # Basic plot
+# Basic plot
 v <- ggplot(faithfuld, aes(waiting, eruptions, z = density))
 v + geom_contour()
 
@@ -187,7 +223,7 @@ ggplot(faithful, aes(waiting, eruptions)) +
 v + geom_contour_filled()
 
 # Setting bins creates evenly spaced contours in the range of the data
-v + geom_contour(bins = 2)
+v + geom_contour(bins = 5)
 v + geom_contour(bins = 10)
 
 # Setting binwidth does the same thing, parameterised by the distance
diff --git a/man/geom_density_2d.Rd b/man/geom_density_2d.Rd
index 86b5f4f855..65623f8260 100644
--- a/man/geom_density_2d.Rd
+++ b/man/geom_density_2d.Rd
@@ -3,16 +3,21 @@
 \name{geom_density_2d}
 \alias{geom_density_2d}
 \alias{geom_density2d}
+\alias{geom_density_2d_filled}
+\alias{geom_density2d_filled}
 \alias{stat_density_2d}
 \alias{stat_density2d}
-\title{Contours of a 2d density estimate}
+\alias{stat_density_2d_filled}
+\alias{stat_density2d_filled}
+\title{Contours of a 2D density estimate}
 \usage{
 geom_density_2d(
   mapping = NULL,
   data = NULL,
-  stat = "density2d",
+  stat = "density_2d",
   position = "identity",
   ...,
+  contour_var = "density",
   lineend = "butt",
   linejoin = "round",
   linemitre = 10,
@@ -21,6 +26,18 @@ geom_density_2d(
   inherit.aes = TRUE
 )
 
+geom_density_2d_filled(
+  mapping = NULL,
+  data = NULL,
+  stat = "density_2d_filled",
+  position = "identity",
+  ...,
+  contour_var = "density",
+  na.rm = FALSE,
+  show.legend = NA,
+  inherit.aes = TRUE
+)
+
 stat_density_2d(
   mapping = NULL,
   data = NULL,
@@ -28,6 +45,23 @@ stat_density_2d(
   position = "identity",
   ...,
   contour = TRUE,
+  contour_var = "density",
+  n = 100,
+  h = NULL,
+  adjust = c(1, 1),
+  na.rm = FALSE,
+  show.legend = NA,
+  inherit.aes = TRUE
+)
+
+stat_density_2d_filled(
+  mapping = NULL,
+  data = NULL,
+  geom = "density_2d_filled",
+  position = "identity",
+  ...,
+  contour = TRUE,
+  contour_var = "density",
   n = 100,
   h = NULL,
   adjust = c(1, 1),
@@ -65,6 +99,10 @@ often aesthetics, used to set an aesthetic to a fixed value, like
 \code{colour = "red"} or \code{size = 3}. They may also be parameters
 to the paired geom/stat.}
 
+\item{contour_var}{Character string identifying the variable to contour
+by. Can be one of \code{"density"}, \code{"ndensity"}, or \code{"count"}. See the section
+on computed variables for details.}
+
 \item{lineend}{Line end style (round, butt, square).}
 
 \item{linejoin}{Line join style (round, mitre, bevel).}
@@ -89,9 +127,9 @@ the default plot specification, e.g. \code{\link[=borders]{borders()}}.}
 \code{geom_density_2d} and \code{stat_density_2d}.}
 
 \item{contour}{If \code{TRUE}, contour the results of the 2d density
-estimation}
+estimation.}
 
-\item{n}{number of grid points in each direction}
+\item{n}{Number of grid points in each direction.}
 
 \item{h}{Bandwidth (vector of length two). If \code{NULL}, estimated
 using \code{\link[MASS:bandwidth.nrd]{MASS::bandwidth.nrd()}}.}
@@ -104,7 +142,9 @@ use half of the default bandwidth.}
 \description{
 Perform a 2D kernel density estimation using \code{\link[MASS:kde2d]{MASS::kde2d()}} and
 display the results with contours. This can be useful for dealing with
-overplotting. This is a 2d version of \code{\link[=geom_density]{geom_density()}}.
+overplotting. This is a 2D version of \code{\link[=geom_density]{geom_density()}}. \code{geom_density_2d()}
+draws contour lines, and \code{geom_density_2d_filled()} draws filled contour
+bands.
 }
 \section{Aesthetics}{
 
@@ -119,17 +159,43 @@ overplotting. This is a 2d version of \code{\link[=geom_density]{geom_density()}
 \item \code{size}
 }
 Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
+
+
+\code{geom_density_2d_filled()} understands the following aesthetics (required aesthetics are in bold):
+\itemize{
+\item \strong{\code{x}}
+\item \strong{\code{y}}
+\item \code{alpha}
+\item \code{colour}
+\item \code{fill}
+\item \code{group}
+\item \code{linetype}
+\item \code{size}
+\item \code{subgroup}
+}
+Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
 }
 
 \section{Computed variables}{
 
-Same as \code{\link[=stat_contour]{stat_contour()}}
-
-With the addition of:
+\code{stat_density_2d()} and \code{stat_density_2d_filled()} compute different
+variables depending on whether contouring is turned on or off. With
+contouring off (\code{contour = FALSE}), both stats behave the same, and the
+following variables are provided:
 \describe{
-\item{density}{the density estimate}
-\item{ndensity}{density estimate, scaled to maximum of 1}
+\item{\code{density}}{The density estimate.}
+\item{\code{ndensity}}{Density estimate, scaled to a maximum of 1.}
+\item{\code{count}}{Density estimate * number of observations in group.}
+\item{\code{n}}{Number of observations in each group.}
 }
+
+With contouring on (\code{contour = TRUE}), either \code{\link[=stat_contour]{stat_contour()}} or
+\code{\link[=stat_contour_filled]{stat_contour_filled()}} (for contour lines or contour bands,
+respectively) is run after the density estimate has been obtained,
+and the computed variables are determined by these stats.
+Contours are calculated for one of the three types of density estimates
+obtained before contouring, \code{density}, \code{ndensity}, and \code{count}. Which
+of those should be used is determined by the \code{contour_var} parameter.
 }
 
 \examples{
@@ -137,9 +203,17 @@ m <- ggplot(faithful, aes(x = eruptions, y = waiting)) +
  geom_point() +
  xlim(0.5, 6) +
  ylim(40, 110)
+
+# contour lines
 m + geom_density_2d()
+
 \donttest{
-m + stat_density_2d(aes(fill = after_stat(level)), geom = "polygon")
+# contour bands
+m + geom_density_2d_filled(alpha = 0.5)
+
+# contour bands and contour lines
+m + geom_density_2d_filled(alpha = 0.5) +
+  geom_density_2d(size = 0.25, colour = "black")
 
 set.seed(4393)
 dsmall <- diamonds[sample(nrow(diamonds), 1000), ]
@@ -148,23 +222,28 @@ d <- ggplot(dsmall, aes(x, y))
 # set of contours for each value of that variable
 d + geom_density_2d(aes(colour = cut))
 
-# Similarly, if you apply faceting to the plot, contours will be
-# drawn for each facet, but the levels will calculated across all facets
-d + stat_density_2d(aes(fill = after_stat(level)), geom = "polygon") +
-  facet_grid(. ~ cut) + scale_fill_viridis_c()
-# To override this behavior (for instace, to better visualize the density
-# within each facet), use after_stat(nlevel)
-d + stat_density_2d(aes(fill = after_stat(nlevel)), geom = "polygon") +
-  facet_grid(. ~ cut) + scale_fill_viridis_c()
-
-# If we turn contouring off, we can use use geoms like tiles:
-d + stat_density_2d(geom = "raster", aes(fill = after_stat(density)), contour = FALSE)
+# If you draw filled contours across multiple facets, the same bins are
+# used across all facets
+d + geom_density_2d_filled() + facet_wrap(vars(cut))
+# If you want to make sure the peak intensity is the same in each facet,
+# use `contour_var = "ndensity"`.
+d + geom_density_2d_filled(contour_var = "ndensity") + facet_wrap(vars(cut))
+# If you want to scale intensity by the number of observations in each group,
+# use `contour_var = "count"`.
+d + geom_density_2d_filled(contour_var = "count") + facet_wrap(vars(cut))
+
+# If we turn contouring off, we can use other geoms, such as tiles:
+d + stat_density_2d(
+  geom = "raster",
+  aes(fill = after_stat(density)),
+  contour = FALSE
+) + scale_fill_viridis_c()
 # Or points:
 d + stat_density_2d(geom = "point", aes(size = after_stat(density)), n = 20, contour = FALSE)
 }
 }
 \seealso{
-\code{\link[=geom_contour]{geom_contour()}} for information about how contours
-are drawn; \code{\link[=geom_bin2d]{geom_bin2d()}} for another way of dealing with
+\code{\link[=geom_contour]{geom_contour()}}, \code{\link[=geom_contour_filled]{geom_contour_filled()}} for information about
+how contours are drawn; \code{\link[=geom_bin2d]{geom_bin2d()}} for another way of dealing with
 overplotting.
 }
diff --git a/man/ggplot2-ggproto.Rd b/man/ggplot2-ggproto.Rd
index 8f4a81ab1b..8fda281c7c 100644
--- a/man/ggplot2-ggproto.Rd
+++ b/man/ggplot2-ggproto.Rd
@@ -58,6 +58,7 @@
 \alias{GeomLine}
 \alias{GeomStep}
 \alias{GeomContour}
+\alias{GeomContourFilled}
 \alias{GeomCrossbar}
 \alias{GeomSegment}
 \alias{GeomCurve}
@@ -65,6 +66,7 @@
 \alias{GeomArea}
 \alias{GeomDensity}
 \alias{GeomDensity2d}
+\alias{GeomDensity2dFilled}
 \alias{GeomDotplot}
 \alias{GeomErrorbar}
 \alias{GeomErrorbarh}
@@ -112,6 +114,7 @@
 \alias{StatContourFilled}
 \alias{StatCount}
 \alias{StatDensity2d}
+\alias{StatDensity2dFilled}
 \alias{StatDensity}
 \alias{StatEcdf}
 \alias{StatEllipse}
diff --git a/tests/testthat/test-stat-density2d.R b/tests/testthat/test-stat-density2d.R
index 340cfa3ccb..3ac5494c90 100644
--- a/tests/testthat/test-stat-density2d.R
+++ b/tests/testthat/test-stat-density2d.R
@@ -16,14 +16,82 @@ test_that("uses scale limits, not data limits", {
   expect_true(max(ret$y) > 35)
 })
 
-# Visual tests --------------------------------------
-
 test_that("stat_density2d can produce contour and raster data", {
   p <- ggplot(faithful, aes(x = eruptions, y = waiting))
 
-  p_contour <- p + stat_density_2d()
+  p_contour_lines <- p + stat_density_2d()
+  p_contour_bands <- p + stat_density_2d_filled()
   p_raster <- p + stat_density_2d(contour = FALSE)
 
-  expect_true("level" %in% names(layer_data(p_contour)))
-  expect_true("density" %in% names(layer_data(p_raster)))
+  d_lines <- layer_data(p_contour_lines)
+  expect_true("level" %in% names(d_lines))
+  expect_false("level_low" %in% names(d_lines))
+  expect_true(is.numeric(d_lines$level))
+
+  d_bands <- layer_data(p_contour_bands)
+  expect_true("level" %in% names(d_bands))
+  expect_true("level_low" %in% names(d_bands))
+  expect_true(is.ordered(d_bands$level))
+
+  d_raster <- layer_data(p_raster)
+  expect_true("density" %in% names(d_raster))
+  expect_true("ndensity" %in% names(d_raster))
+  expect_true("count" %in% names(d_raster))
+  expect_true(unique(d_raster$level) == 1)
+  expect_true(unique(d_raster$piece) == 1)
+
+  # stat_density_2d() and stat_density_2d_filled() produce identical
+  # density output with `contour = FALSE`
+  # (`fill` and `colour` will differ due to different default aesthetic mappings)
+  d_raster2 <- layer_data(p + stat_density_2d_filled(contour = FALSE))
+  expect_identical(d_raster$x, d_raster2$x)
+  expect_identical(d_raster$y, d_raster2$y)
+  expect_identical(d_raster$density, d_raster2$density)
+  expect_identical(d_raster$ndensity, d_raster2$ndensity)
+  expect_identical(d_raster$count, d_raster2$count)
+
+  # stat_density_2d() with contouring is the same as stat_contour() on calculated density
+  p_lines2 <- ggplot(d_raster, aes(x, y, z = density)) + stat_contour()
+  d_lines2 <- layer_data(p_lines2)
+  expect_identical(d_lines$x, d_lines2$x)
+  expect_identical(d_lines$y, d_lines2$y)
+  expect_identical(d_lines$piece, d_lines2$piece)
+  expect_identical(d_lines$group, d_lines2$group)
+  expect_identical(d_lines$level, d_lines2$level)
+
+  # same for stat_density_2d_filled()
+  p_bands2 <- ggplot(d_raster, aes(x, y, z = density)) + stat_contour_filled()
+  d_bands2 <- layer_data(p_bands2)
+  expect_identical(d_bands$x, d_bands2$x)
+  expect_identical(d_bands$y, d_bands2$y)
+  expect_identical(d_bands$piece, d_bands2$piece)
+  expect_identical(d_bands$group, d_bands2$group)
+  expect_identical(d_bands$level, d_bands2$level)
+  expect_identical(d_bands$level_mid, d_bands2$level_mid)
+
+  # and for contour_var = "ndensity"
+  p_contour_lines <- p + stat_density_2d(contour_var = "ndensity")
+  d_lines <- layer_data(p_contour_lines)
+  p_lines2 <- ggplot(d_raster, aes(x, y, z = ndensity)) + stat_contour()
+  d_lines2 <- layer_data(p_lines2)
+  expect_identical(d_lines$x, d_lines2$x)
+  expect_identical(d_lines$y, d_lines2$y)
+  expect_identical(d_lines$piece, d_lines2$piece)
+  expect_identical(d_lines$group, d_lines2$group)
+  expect_identical(d_lines$level, d_lines2$level)
+
+  # and for contour_var = "count"
+  p_contour_bands <- p + stat_density_2d_filled(contour_var = "count")
+  d_bands <- layer_data(p_contour_bands)
+  p_bands2 <- ggplot(d_raster, aes(x, y, z = count)) + stat_contour_filled()
+  d_bands2 <- layer_data(p_bands2)
+  expect_identical(d_bands$x, d_bands2$x)
+  expect_identical(d_bands$y, d_bands2$y)
+  expect_identical(d_bands$piece, d_bands2$piece)
+  expect_identical(d_bands$group, d_bands2$group)
+  expect_identical(d_bands$level, d_bands2$level)
+  expect_identical(d_bands$level_mid, d_bands2$level_mid)
+
+  # error on incorrect contouring variable
+  expect_error(ggplot_build(p + stat_density_2d(contour_var = "abcd")))
 })