diff --git a/DESCRIPTION b/DESCRIPTION
index 20143da6b5..f07e81ac0e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -250,9 +250,10 @@ Collate:
     'utilities-resolution.r'
     'utilities-table.r'
     'utilities-tidy-eval.R'
+    'with-sf.R'
     'zxx.r'
     'zzz.r'
 VignetteBuilder: knitr
-RoxygenNote: 7.0.1
+RoxygenNote: 7.0.2
 Roxygen: list(markdown = TRUE)
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index af5d41396e..40c1fffe9f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -630,6 +630,7 @@ export(update_labels)
 export(update_stat_defaults)
 export(vars)
 export(waiver)
+export(with_sf)
 export(wrap_dims)
 export(xlab)
 export(xlim)
diff --git a/R/with-sf.R b/R/with-sf.R
new file mode 100644
index 0000000000..4439828000
--- /dev/null
+++ b/R/with-sf.R
@@ -0,0 +1,116 @@
+#' Turn regular geoms/stats into sf-ified geoms/stats
+#'
+#' This function allows you to wrap ordinary geoms and stats for use within
+#' a `geom_sf()`/`coord_sf()` framework. By default, it assumes that coordinates
+#' in the wrapped geoms/stats are given in standard latitude/longitude coordinates,
+#' and it transforms those coordinates into those used by `coord_sf()`. However,
+#' it is also possible to specify an alternative coordinate reference system (crs)
+#' for the geoms/stats that are being wrapped.
+#' @examples
+#' \donttest{
+#' nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
+#'
+#' df <- data.frame(lat = 35.76667, long = -78.63333, name = "Raleigh")
+#'
+#' ggplot(df, aes(x = long, y = lat)) +
+#'   geom_sf(data = nc, size = 0.1, fill = "white", inherit.aes = FALSE) +
+#'   with_sf(
+#'     geom_point(),
+#'     geom_text(aes(label = name), hjust = -.1, vjust = 1.1)
+#'   ) +
+#'   coord_sf(crs = 2264)
+#'
+#' ggplot(nc) +
+#'   geom_sf(size = 0.1, fill = "white") +
+#'   with_sf(
+#'     geom_hline(yintercept = 34:36),
+#'     geom_vline(xintercept = c(-84, -80, -76))
+#'   ) +
+#'   coord_sf(crs = 2264)
+#' }
+#' @param ... One or more geom/stat or lists of geoms/stats to be transformed
+#' @param crs Coordinate reference system of the origin layer. Defaults to
+#'   longitude and latitude in WGS84 (EPSG:4326), the World Geodetic System
+#'   from 1984.
+#' @export
+with_sf <- function(..., crs = 4326) {
+  x <- list(...)
+  with_sf_impl(x, crs)
+}
+
+with_sf_impl <- function(x, crs, ...) {
+  UseMethod("with_sf_impl")
+}
+
+with_sf_impl.default <- function(x, crs, ...) {
+  stop(
+    "Can't convert object of class `", class(x), "` to an sf layer.\n",
+    call. = FALSE
+  )
+}
+
+with_sf_impl.list <- function(x, crs, ...) {
+  l <- list()
+  # for some reason lapply() version of this doesn't work
+  for (i in seq_along(x)) {
+    l[[i]] <- with_sf_impl(x[[i]], crs)
+  }
+  if (length(l) == 1) return(l[[1]])
+  l
+}
+
+with_sf_impl.Layer <- function(x, crs, ...) {
+  parent_geom <- x$geom
+  ggproto(NULL, x,
+    geom = ggproto('SfifiedGeom', parent_geom,
+      draw_panel = function(data, panel_params, coord, ...) {
+        sfified_coord <- ggproto('SfifiedCoord', coord,
+          transform = function(data, range) {
+            data <- sfified_transform(data, crs = coord$crs, source_crs = crs)
+            coord$transform(data, range)
+          },
+
+          backtransform_range = function(panel_params) {
+            # we create a bounding box, transform all four corners, and then extract the max extent
+            x <- panel_params$x_range
+            y <- panel_params$y_range
+            data <- list(x = c(x, x), y = c(y, rev(y)))
+            data <- sfified_transform(data, crs = crs, source_crs = coord$crs)
+            list(x = range(data$x), y = range(data$y))
+          },
+
+          # normally coord_sf() pretends to be linear, but as used here it is not
+          is_linear = function() FALSE,
+
+          distance = function(self, x, y, panel_params) {
+            d <- self$backtransform_range(panel_params)
+            max_dist <- dist_euclidean(d$x, d$y)
+            dist_euclidean(x, y) / max_dist
+          }
+        )
+        parent_geom$draw_panel(data, panel_params, sfified_coord, ...)
+      },
+
+      parameters = function(self, extra = FALSE) {
+        # make sure we extract parameters of wrapped geom correctly
+        parent_geom$parameters(extra)
+      }
+    )
+  )
+}
+
+
+sfified_transform <- function(data, crs = NULL, source_crs = 4326) {
+  if (is.null(crs)) {
+    return(data)
+  }
+
+  sf_data <- sf::st_sfc(
+    sf::st_multipoint(cbind(data$x, data$y)),
+    crs = source_crs
+  )
+  sf_data_trans <- sf::st_transform(sf_data, crs)[[1]]
+  data$x <- sf_data_trans[, 1]
+  data$y <- sf_data_trans[, 2]
+  data
+}
diff --git a/man/with_sf.Rd b/man/with_sf.Rd
new file mode 100644
index 0000000000..eb8dfee017
--- /dev/null
+++ b/man/with_sf.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/with-sf.R
+\name{with_sf}
+\alias{with_sf}
+\title{Turn regular geoms/stats into sf-ified geoms/stats}
+\usage{
+with_sf(..., crs = 4326)
+}
+\arguments{
+\item{...}{One or more geom/stat or lists of geoms/stats to be transformed}
+
+\item{crs}{Coordinate reference system of the origin layer. Defaults to
+longitude and latitude in WGS84 (EPSG:4326), the World Geodetic System
+from 1984.}
+}
+\description{
+This function allows you to wrap ordinary geoms and stats for use within
+a \code{geom_sf()}/\code{coord_sf()} framework. By default, it assumes that coordinates
+in the wrapped geoms/stats are given in standard latitude/longitude coordinates,
+and it transforms those coordinates into those used by \code{coord_sf()}. However,
+it is also possible to specify an alternative coordinate reference system (crs)
+for the geoms/stats that are being wrapped.
+}
+\examples{
+\donttest{
+nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
+
+df <- data.frame(lat = 35.76667, long = -78.63333, name = "Raleigh")
+
+ggplot(df, aes(x = long, y = lat)) +
+  geom_sf(data = nc, size = 0.1, fill = "white", inherit.aes = FALSE) +
+  with_sf(
+    geom_point(),
+    geom_text(aes(label = name), hjust = -.1, vjust = 1.1)
+  ) +
+  coord_sf(crs = 2264)
+
+ggplot(nc) +
+  geom_sf(size = 0.1, fill = "white") +
+  with_sf(
+    geom_hline(yintercept = 34:36),
+    geom_vline(xintercept = c(-84, -80, -76))
+  ) +
+  coord_sf(crs = 2264)
+}
+}