diff --git a/R/coord-map.r b/R/coord-map.r
index 4f55880f63..d8f4624186 100644
--- a/R/coord-map.r
+++ b/R/coord-map.r
@@ -1,5 +1,8 @@
 #' Map projections
 #'
+#' @description
+#' `r lifecycle::badge("superseded")`
+#'
 #' `coord_map()` projects a portion of the earth, which is approximately
 #' spherical, onto a flat 2D plane using any projection defined by the
 #' `mapproj` package. Map projections do not, in general, preserve straight
@@ -7,7 +10,15 @@
 #' quick approximation that does preserve straight lines. It works best for
 #' smaller areas closer to the equator.
 #'
-#' In general, map projections must account for the fact that the actual length
+#' Both `coord_map()` and `coord_quickmap()`
+#' are superseded by [`coord_sf()`], and should no longer be used in new
+#' code. All regular (non-sf) geoms can be used with `coord_sf()` by
+#' setting the default coordinate system via the `default_crs` argument.
+#' See also the examples for [`annotation_map()`] and [`geom_map()`].
+#'
+#' @details
+#'
+#' Map projections must account for the fact that the actual length
 #' (in km) of one degree of longitude varies between the equator and the pole.
 #' Near the equator, the ratio between the lengths of one degree of latitude and
 #' one degree of longitude is approximately 1. Near the pole, it tends
diff --git a/R/geom-map.r b/R/geom-map.r
index c04573e2eb..c26178e4ac 100644
--- a/R/geom-map.r
+++ b/R/geom-map.r
@@ -3,7 +3,11 @@ NULL
 
 #' Polygons from a reference map
 #'
-#' This is pure annotation, so does not affect position scales.
+#' Display polygons as a map. This is meant as annotation, so it does not
+#' affect position scales. Note that this function predates the [`geom_sf()`]
+#' framework and does not work with sf geometry columns as input. However,
+#' it can be used in conjunction with `geom_sf()` layers and/or
+#' [`coord_sf()`] (see examples).
 #'
 #' @eval rd_aesthetics("geom", "map")
 #' @export
@@ -14,10 +18,12 @@ NULL
 #' @inheritParams layer
 #' @inheritParams geom_point
 #' @examples
-#' # When using geom_polygon, you will typically need two data frames:
-#' # one contains the coordinates of each polygon (positions),  and the
-#' # other the values associated with each polygon (values).  An id
-#' # variable links the two together
+#' # First, a made-up example containing a few polygons, to explain
+#' # how `geom_map()` works. It requires two data frames:
+#' # One contains the coordinates of each polygon (`positions`), and is
+#' # provided via the `map` argument. The other contains the
+#' # other the values associated with each polygon (`values`).  An id
+#' # variable links the two together.
 #'
 #' ids <- factor(c("1.1", "2.1", "1.2", "2.2", "1.3", "2.3"))
 #'
@@ -44,7 +50,7 @@ NULL
 #'   geom_map(aes(map_id = id), map = positions) +
 #'   expand_limits(positions) + ylim(0, 3)
 #'
-#' # Better example
+#' # Now some examples with real maps
 #' if (require(maps)) {
 #'
 #'   crimes <- data.frame(state = tolower(rownames(USArrests)), USArrests)
@@ -56,15 +62,32 @@ NULL
 #'   crimes_long <- do.call("rbind", vars)
 #'
 #'   states_map <- map_data("state")
+#'
+#'   # without geospatial coordinate system, the resulting plot
+#'   # looks weird
 #'   ggplot(crimes, aes(map_id = state)) +
 #'     geom_map(aes(fill = Murder), map = states_map) +
 #'     expand_limits(x = states_map$long, y = states_map$lat)
 #'
-#'   last_plot() + coord_map()
-#'   ggplot(crimes_long, aes(map_id = state)) +
-#'     geom_map(aes(fill = value), map = states_map) +
-#'     expand_limits(x = states_map$long, y = states_map$lat) +
-#'     facet_wrap( ~ variable)
+#'   # in combination with `coord_sf()` we get an appropriate result
+#'   ggplot(crimes, aes(map_id = state)) +
+#'     geom_map(aes(fill = Murder), map = states_map) +
+#'     # crs = 5070 is a Conus Albers projection for North America,
+#'     #   see: https://epsg.io/5070
+#'     # default_crs = 4326 tells coord_sf() that the input map data
+#'     #   are in longitude-latitude format
+#'     coord_sf(
+#'       crs = 5070, default_crs = 4326,
+#'       xlim = c(-125, -70), ylim = c(25, 52)
+#'     )
+#'
+#'  ggplot(crimes_long, aes(map_id = state)) +
+#'    geom_map(aes(fill = value), map = states_map) +
+#'    coord_sf(
+#'      crs = 5070, default_crs = 4326,
+#'      xlim = c(-125, -70), ylim = c(25, 52)
+#'    ) +
+#'    facet_wrap(~variable)
 #' }
 geom_map <- function(mapping = NULL, data = NULL,
                      stat = "identity",
diff --git a/man/coord_map.Rd b/man/coord_map.Rd
index c3ead2cdf6..8fb534af48 100644
--- a/man/coord_map.Rd
+++ b/man/coord_map.Rd
@@ -43,15 +43,23 @@ the limits to ensure that data and axes don't overlap. If \code{FALSE},
 limits are taken exactly from the data or \code{xlim}/\code{ylim}.}
 }
 \description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#superseded}{\figure{lifecycle-superseded.svg}{options: alt='[Superseded]'}}}{\strong{[Superseded]}}
+
 \code{coord_map()} projects a portion of the earth, which is approximately
 spherical, onto a flat 2D plane using any projection defined by the
 \code{mapproj} package. Map projections do not, in general, preserve straight
 lines, so this requires considerable computation. \code{coord_quickmap()} is a
 quick approximation that does preserve straight lines. It works best for
 smaller areas closer to the equator.
+
+Both \code{coord_map()} and \code{coord_quickmap()}
+are superseded by \code{\link[=coord_sf]{coord_sf()}}, and should no longer be used in new
+code. All regular (non-sf) geoms can be used with \code{coord_sf()} by
+setting the default coordinate system via the \code{default_crs} argument.
+See also the examples for \code{\link[=annotation_map]{annotation_map()}} and \code{\link[=geom_map]{geom_map()}}.
 }
 \details{
-In general, map projections must account for the fact that the actual length
+Map projections must account for the fact that the actual length
 (in km) of one degree of longitude varies between the equator and the pole.
 Near the equator, the ratio between the lengths of one degree of latitude and
 one degree of longitude is approximately 1. Near the pole, it tends
diff --git a/man/geom_map.Rd b/man/geom_map.Rd
index a992619015..18e765ba8e 100644
--- a/man/geom_map.Rd
+++ b/man/geom_map.Rd
@@ -66,7 +66,11 @@ that define both data and aesthetics and shouldn't inherit behaviour from
 the default plot specification, e.g. \code{\link[=borders]{borders()}}.}
 }
 \description{
-This is pure annotation, so does not affect position scales.
+Display polygons as a map. This is meant as annotation, so it does not
+affect position scales. Note that this function predates the \code{\link[=geom_sf]{geom_sf()}}
+framework and does not work with sf geometry columns as input. However,
+it can be used in conjunction with \code{geom_sf()} layers and/or
+\code{\link[=coord_sf]{coord_sf()}} (see examples).
 }
 \section{Aesthetics}{
 
@@ -85,10 +89,12 @@ Learn more about setting these aesthetics in \code{vignette("ggplot2-specs")}.
 }
 
 \examples{
-# When using geom_polygon, you will typically need two data frames:
-# one contains the coordinates of each polygon (positions),  and the
-# other the values associated with each polygon (values).  An id
-# variable links the two together
+# First, a made-up example containing a few polygons, to explain
+# how `geom_map()` works. It requires two data frames:
+# One contains the coordinates of each polygon (`positions`), and is
+# provided via the `map` argument. The other contains the
+# other the values associated with each polygon (`values`).  An id
+# variable links the two together.
 
 ids <- factor(c("1.1", "2.1", "1.2", "2.2", "1.3", "2.3"))
 
@@ -115,7 +121,7 @@ ggplot(values, aes(fill = value)) +
   geom_map(aes(map_id = id), map = positions) +
   expand_limits(positions) + ylim(0, 3)
 
-# Better example
+# Now some examples with real maps
 if (require(maps)) {
 
   crimes <- data.frame(state = tolower(rownames(USArrests)), USArrests)
@@ -127,14 +133,31 @@ if (require(maps)) {
   crimes_long <- do.call("rbind", vars)
 
   states_map <- map_data("state")
+
+  # without geospatial coordinate system, the resulting plot
+  # looks weird
   ggplot(crimes, aes(map_id = state)) +
     geom_map(aes(fill = Murder), map = states_map) +
     expand_limits(x = states_map$long, y = states_map$lat)
 
-  last_plot() + coord_map()
-  ggplot(crimes_long, aes(map_id = state)) +
-    geom_map(aes(fill = value), map = states_map) +
-    expand_limits(x = states_map$long, y = states_map$lat) +
-    facet_wrap( ~ variable)
+  # in combination with `coord_sf()` we get an appropriate result
+  ggplot(crimes, aes(map_id = state)) +
+    geom_map(aes(fill = Murder), map = states_map) +
+    # crs = 5070 is a Conus Albers projection for North America,
+    #   see: https://epsg.io/5070
+    # default_crs = 4326 tells coord_sf() that the input map data
+    #   are in longitude-latitude format
+    coord_sf(
+      crs = 5070, default_crs = 4326,
+      xlim = c(-125, -70), ylim = c(25, 52)
+     )
+
+ ggplot(crimes_long, aes(map_id = state)) +
+   geom_map(aes(fill = value), map = states_map) +
+   coord_sf(
+     crs = 5070, default_crs = 4326,
+     xlim = c(-125, -70), ylim = c(25, 52)
+   ) +
+   facet_wrap(~variable)
 }
 }