diff --git a/DESCRIPTION b/DESCRIPTION
index 349f905e59..3653d3a6ed 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -171,6 +171,7 @@ Collate:
     'geom-violin.R'
     'geom-vline.R'
     'ggplot2-package.R'
+    'ggunit.R'
     'grob-absolute.R'
     'grob-dotstack.R'
     'grob-null.R'
@@ -280,5 +281,6 @@ Collate:
     'utilities-patterns.R'
     'utilities-resolution.R'
     'utilities-tidy-eval.R'
+    'utilities-unit.R'
     'zxx.R'
     'zzz.R'
diff --git a/NAMESPACE b/NAMESPACE
index 967573b174..04c7faa034 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,18 +5,24 @@ S3method("$",ggproto_parent)
 S3method("$",theme)
 S3method("$<-",uneval)
 S3method("+",gg)
+S3method("[",ggunit)
 S3method("[",mapped_discrete)
 S3method("[",uneval)
+S3method("[<-",ggunit)
 S3method("[<-",mapped_discrete)
 S3method("[<-",uneval)
 S3method("[[",ggproto)
+S3method("[[<-",ggunit)
 S3method("[[<-",uneval)
 S3method(.DollarNames,ggproto)
+S3method(Ops,ggunit)
+S3method(Summary,ggunit)
 S3method(as.data.frame,mapped_discrete)
 S3method(as.list,ggproto)
 S3method(autolayer,default)
 S3method(autoplot,default)
 S3method(c,mapped_discrete)
+S3method(chooseOpsMethod,ggunit)
 S3method(drawDetails,zeroGrob)
 S3method(element_grob,element_blank)
 S3method(element_grob,element_line)
@@ -85,6 +91,9 @@ S3method(heightDetails,titleGrob)
 S3method(heightDetails,zeroGrob)
 S3method(interleave,default)
 S3method(interleave,unit)
+S3method(is.finite,ggunit)
+S3method(is.infinite,ggunit)
+S3method(is.na,ggunit)
 S3method(limits,Date)
 S3method(limits,POSIXct)
 S3method(limits,POSIXlt)
@@ -92,6 +101,7 @@ S3method(limits,character)
 S3method(limits,factor)
 S3method(limits,numeric)
 S3method(makeContext,dotstackGrob)
+S3method(makeContext,ggplot2_pointsGrob)
 S3method(merge_element,default)
 S3method(merge_element,element)
 S3method(merge_element,element_blank)
@@ -109,9 +119,11 @@ S3method(print,ggplot)
 S3method(print,ggplot2_bins)
 S3method(print,ggproto)
 S3method(print,ggproto_method)
+S3method(print,ggunit)
 S3method(print,rel)
 S3method(print,theme)
 S3method(print,uneval)
+S3method(rescale,unit)
 S3method(scale_type,AsIs)
 S3method(scale_type,Date)
 S3method(scale_type,POSIXt)
@@ -126,27 +138,63 @@ S3method(scale_type,logical)
 S3method(scale_type,numeric)
 S3method(scale_type,ordered)
 S3method(scale_type,sfc)
+S3method(scale_type,unit)
+S3method(sign,ggunit)
 S3method(single_value,default)
 S3method(single_value,factor)
 S3method(summary,ggplot)
 S3method(vec_cast,character.mapped_discrete)
+S3method(vec_cast,double.ggunit)
 S3method(vec_cast,double.mapped_discrete)
 S3method(vec_cast,factor.mapped_discrete)
+S3method(vec_cast,ggunit.character)
+S3method(vec_cast,ggunit.double)
+S3method(vec_cast,ggunit.ggunit)
+S3method(vec_cast,ggunit.integer)
+S3method(vec_cast,ggunit.list)
+S3method(vec_cast,ggunit.logical)
+S3method(vec_cast,ggunit.mapped_discrete)
+S3method(vec_cast,ggunit.simpleUnit)
+S3method(vec_cast,ggunit.unit)
+S3method(vec_cast,integer.ggunit)
 S3method(vec_cast,integer.mapped_discrete)
+S3method(vec_cast,logical.ggunit)
 S3method(vec_cast,mapped_discrete.double)
 S3method(vec_cast,mapped_discrete.factor)
+S3method(vec_cast,mapped_discrete.ggunit)
 S3method(vec_cast,mapped_discrete.integer)
 S3method(vec_cast,mapped_discrete.logical)
 S3method(vec_cast,mapped_discrete.mapped_discrete)
+S3method(vec_cast,simpleUnit.ggunit)
+S3method(vec_cast,simpleUnit.unit)
+S3method(vec_cast,unit.ggunit)
+S3method(vec_cast,unit.simpleUnit)
+S3method(vec_cast,unit.unit)
+S3method(vec_proxy,simpleUnit)
+S3method(vec_proxy,unit)
 S3method(vec_ptype2,character.mapped_discrete)
+S3method(vec_ptype2,double.ggunit)
 S3method(vec_ptype2,double.mapped_discrete)
 S3method(vec_ptype2,factor.mapped_discrete)
+S3method(vec_ptype2,ggunit.double)
+S3method(vec_ptype2,ggunit.ggunit)
+S3method(vec_ptype2,ggunit.integer)
+S3method(vec_ptype2,ggunit.simpleUnit)
+S3method(vec_ptype2,ggunit.unit)
+S3method(vec_ptype2,integer.ggunit)
 S3method(vec_ptype2,integer.mapped_discrete)
 S3method(vec_ptype2,mapped_discrete.character)
 S3method(vec_ptype2,mapped_discrete.double)
 S3method(vec_ptype2,mapped_discrete.factor)
 S3method(vec_ptype2,mapped_discrete.integer)
 S3method(vec_ptype2,mapped_discrete.mapped_discrete)
+S3method(vec_ptype2,simpleUnit.ggunit)
+S3method(vec_ptype2,simpleUnit.unit)
+S3method(vec_ptype2,unit.ggunit)
+S3method(vec_ptype2,unit.simpleUnit)
+S3method(vec_ptype2,unit.unit)
+S3method(vec_restore,ggunit)
+S3method(vec_restore,unit)
 S3method(widthDetails,titleGrob)
 S3method(widthDetails,zeroGrob)
 export("%+%")
diff --git a/R/coord-.R b/R/coord-.R
index 8c4313baf7..26c3f9449b 100644
--- a/R/coord-.R
+++ b/R/coord-.R
@@ -170,7 +170,11 @@ Coord <- ggproto("Coord",
     panel_params
   },
 
-  transform = function(data, range) NULL,
+  transform = function(self, data, range) {
+    .expose_units(self$transform_native(.ignore_units(data), range))
+  },
+
+  transform_native = function(data, range) NULL,
 
   distance = function(x, y, panel_params) NULL,
 
diff --git a/R/coord-cartesian-.R b/R/coord-cartesian-.R
index 74f46433db..5977d4e45a 100644
--- a/R/coord-cartesian-.R
+++ b/R/coord-cartesian-.R
@@ -93,7 +93,7 @@ CoordCartesian <- ggproto("CoordCartesian", Coord,
     self$range(panel_params)
   },
 
-  transform = function(data, panel_params) {
+  transform_native = function(data, panel_params) {
     data <- transform_position(data, panel_params$x$rescale, panel_params$y$rescale)
     transform_position(data, squish_infinite, squish_infinite)
   },
diff --git a/R/coord-flip.R b/R/coord-flip.R
index 1f3848fb8a..c688be7cf0 100644
--- a/R/coord-flip.R
+++ b/R/coord-flip.R
@@ -59,7 +59,7 @@ coord_flip <- function(xlim = NULL, ylim = NULL, expand = TRUE, clip = "on") {
 #' @export
 CoordFlip <- ggproto("CoordFlip", CoordCartesian,
 
-  transform = function(data, panel_params) {
+  transform_native = function(data, panel_params) {
     data <- flip_axis_labels(data)
     CoordCartesian$transform(data, panel_params)
   },
diff --git a/R/coord-map.R b/R/coord-map.R
index ee0f6ad139..649d8f1436 100644
--- a/R/coord-map.R
+++ b/R/coord-map.R
@@ -153,7 +153,7 @@ coord_map <- function(projection="mercator", ..., parameters = NULL, orientation
 #' @export
 CoordMap <- ggproto("CoordMap", Coord,
 
-  transform = function(self, data, panel_params) {
+  transform_native = function(self, data, panel_params) {
     trans <- mproject(self, data$x, data$y, panel_params$orientation)
     out <- cunion(trans[c("x", "y")], data)
 
diff --git a/R/coord-munch.R b/R/coord-munch.R
index 6f2bbb2afb..6f877fc506 100644
--- a/R/coord-munch.R
+++ b/R/coord-munch.R
@@ -15,6 +15,8 @@
 coord_munch <- function(coord, data, range, segment_length = 0.01, is_closed = FALSE) {
   if (coord$is_linear()) return(coord$transform(data, range))
 
+  data <- .ignore_units(data)
+
   if (is_closed) {
     data <- close_poly(data)
   }
@@ -44,7 +46,7 @@ coord_munch <- function(coord, data, range, segment_length = 0.01, is_closed = F
     runs <- vec_run_sizes(munched[, group_cols, drop = FALSE])
     munched <- vec_slice(munched, -(cumsum(runs)))
   }
-  coord$transform(munched, range)
+  coord$transform(.expose_units(munched), range)
 }
 
 # For munching, only grobs are lines and polygons: everything else is
diff --git a/R/coord-polar.R b/R/coord-polar.R
old mode 100644
new mode 100755
index 1e30adcd2b..9d59fa843d
--- a/R/coord-polar.R
+++ b/R/coord-polar.R
@@ -167,7 +167,7 @@ CoordPolar <- ggproto("CoordPolar", Coord,
     panel_params
   },
 
-  transform = function(self, data, panel_params) {
+  transform_native = function(self, data, panel_params) {
     arc  <- self$start + c(0, 2 * pi)
     dir  <- self$direction
     data <- rename_data(self, data)
@@ -316,15 +316,15 @@ rename_data <- function(coord, data) {
 }
 
 theta_rescale_no_clip <- function(x, range, arc = c(0, 2 * pi), direction = 1) {
-  rescale(x, to = arc, from = range) * direction
+  transform_native_units(x, function(x) rescale(x, to = arc, from = range) * direction)
 }
 
 theta_rescale <- function(x, range, arc = c(0, 2 * pi), direction = 1) {
   x <- squish_infinite(x, range)
-  rescale(x, to = arc, from = range) %% (2 * pi) * direction
+  transform_native_units(x, function(x) rescale(x, to = arc, from = range) %% (2 * pi) * direction)
 }
 
 r_rescale <- function(x, range, donut = c(0, 0.4)) {
   x <- squish_infinite(x, range)
-  rescale(x, donut, range)
+  transform_native_units(x, function(x) rescale(x, donut, range))
 }
diff --git a/R/coord-radial.R b/R/coord-radial.R
index 70aa211898..607f94084a 100644
--- a/R/coord-radial.R
+++ b/R/coord-radial.R
@@ -219,7 +219,7 @@ CoordRadial <- ggproto("CoordRadial", Coord,
     panel_params
   },
 
-  transform = function(self, data, panel_params) {
+  transform_native = function(self, data, panel_params) {
     data <- rename_data(self, data)
     bbox <- panel_params$bbox %||% list(x = c(0, 1), y = c(0, 1))
     arc  <- panel_params$arc  %||% c(0, 2 * pi)
diff --git a/R/coord-sf.R b/R/coord-sf.R
index 331ca4f1f0..71a666e326 100644
--- a/R/coord-sf.R
+++ b/R/coord-sf.R
@@ -78,7 +78,7 @@ CoordSf <- ggproto("CoordSf", CoordCartesian,
     self$params$bbox <- bbox
   },
 
-  transform = function(self, data, panel_params) {
+  transform_native = function(self, data, panel_params) {
     # we need to transform all non-sf data into the correct coordinate system
     source_crs <- panel_params$default_crs
     target_crs <- panel_params$crs
diff --git a/R/coord-transform.R b/R/coord-transform.R
index 79d651e8af..3fc8f36816 100644
--- a/R/coord-transform.R
+++ b/R/coord-transform.R
@@ -127,7 +127,7 @@ CoordTrans <- ggproto("CoordTrans", Coord,
     )
   },
 
-  transform = function(self, data, panel_params) {
+  transform_native = function(self, data, panel_params) {
     # trans_x() and trans_y() needs to keep Inf values because this can be called
     # in guide_transform.axis()
     trans_x <- function(data) {
diff --git a/R/fortify-map.R b/R/fortify-map.R
index d0dc76b716..6237ffebee 100644
--- a/R/fortify-map.R
+++ b/R/fortify-map.R
@@ -37,7 +37,7 @@ fortify.map <- function(model, data, ...) {
   names <- inject(rbind(!!!names))
   df$region <- names[df$group, 1]
   df$subregion <- names[df$group, 2]
-  df[stats::complete.cases(df$lat, df$long), ]
+  df[vec_detect_complete(df$lat, df$long), ]
 }
 
 #' Create a data frame of map data
diff --git a/R/geom-path.R b/R/geom-path.R
index cf9e59976c..696f680031 100644
--- a/R/geom-path.R
+++ b/R/geom-path.R
@@ -138,7 +138,7 @@ GeomPath <- ggproto("GeomPath", Geom,
     # Drop missing values at the start or end of a line - can't drop in the
     # middle since you expect those to be shown by a break in the line
     aesthetics <- c(self$required_aes, self$non_missing_aes)
-    complete <- stats::complete.cases(data[names(data) %in% aesthetics])
+    complete <- vec_detect_complete(data[names(data) %in% aesthetics])
     kept <- stats::ave(complete, data$group, FUN = keep_mid_true)
     data <- data[kept, ]
 
diff --git a/R/geom-point.R b/R/geom-point.R
old mode 100644
new mode 100755
index 1b39a11d46..f014853eac
--- a/R/geom-point.R
+++ b/R/geom-point.R
@@ -121,21 +121,26 @@ GeomPoint <- ggproto("GeomPoint", Geom,
     if (is.character(data$shape)) {
       data$shape <- translate_shape_string(data$shape)
     }
-
     coords <- coord$transform(data, panel_params)
+
     stroke_size <- coords$stroke
-    stroke_size[is.na(stroke_size)] <- 0
+    if (!is.unit(stroke_size)) stroke_size <- unit(stroke_size * .stroke, "pt")
+    native_units(stroke_size) <- rescale(native_units(stroke_size), from = c(0, diff(coord$range(panel_params)$x)))
+    stroke_size[is.na(stroke_size)] <- unit(0, "pt")
+
+    font_size <- coords$size
+    if (!is.unit(font_size)) font_size <- unit(font_size * .pt, "pt")
+    native_units(font_size) <- rescale(native_units(font_size), from = c(0, diff(coord$range(panel_params)$x)))
+
     ggname("geom_point",
-      pointsGrob(
+      ggplot2_pointsGrob(
         coords$x, coords$y,
         pch = coords$shape,
-        gp = gpar(
-          col = alpha(coords$colour, coords$alpha),
-          fill = fill_alpha(coords$fill, coords$alpha),
-          # Stroke is added around the outside of the point
-          fontsize = coords$size * .pt + stroke_size * .stroke / 2,
-          lwd = coords$stroke * .stroke / 2
-        )
+        col = alpha(coords$colour, coords$alpha),
+        fill = fill_alpha(coords$fill, coords$alpha),
+        # Stroke is added around the outside of the point
+        fontsize = font_size + stroke_size / 2,
+        lwd = stroke_size / 2
       )
     )
   },
@@ -143,6 +148,33 @@ GeomPoint <- ggproto("GeomPoint", Geom,
   draw_key = draw_key_point
 )
 
+ggplot2_pointsGrob <- function(
+  x, y, pch = 1, vp = NULL,
+  fontsize = 12, lwd = 1, col = "black", fill = "white"
+) {
+  grob(
+    x = x, y = y, pch = pch, vp = vp,
+    fontsize = fontsize, lwd = lwd, col = col, fill = fill,
+    cl = "ggplot2_pointsGrob"
+  )
+}
+
+#' @export
+makeContext.ggplot2_pointsGrob <- function(x) {
+  pointsGrob(
+    x$x,
+    x$y,
+    pch = x$pch,
+    gp = gpar(
+      col = x$col,
+      fill = x$fill,
+      # Stroke is added around the outside of the point
+      fontsize = convertUnit(x$fontsize + x$lwd, unitTo = "pt", valueOnly = TRUE),
+      lwd = convertUnit(x$lwd, unitTo = "pt", valueOnly = TRUE)
+    )
+  )
+}
+
 #' Translating shape strings
 #'
 #' `translate_shape_string()` is a helper function for translating point shapes
diff --git a/R/geom-ribbon.R b/R/geom-ribbon.R
index d93df77850..18569f4335 100644
--- a/R/geom-ribbon.R
+++ b/R/geom-ribbon.R
@@ -131,7 +131,7 @@ GeomRibbon <- ggproto("GeomRibbon", Geom,
                         flipped_aes = FALSE, outline.type = "both") {
     data <- check_linewidth(data, snake_class(self))
     data <- flip_data(data, flipped_aes)
-    if (na.rm) data <- data[stats::complete.cases(data[c("x", "ymin", "ymax")]), ]
+    if (na.rm) data <- data[vec_detect_complete(data[c("x", "ymin", "ymax")]), ]
     data <- data[order(data$group), ]
 
     # Check that aesthetics are constant
@@ -148,7 +148,7 @@ GeomRibbon <- ggproto("GeomRibbon", Geom,
     # has distinct polygon numbers for sequences of non-NA values and NA
     # for NA values in the original data.  Example: c(NA, 2, 2, 2, NA, NA,
     # 4, 4, 4, NA)
-    missing_pos <- !stats::complete.cases(data[c("x", "ymin", "ymax")])
+    missing_pos <- !vec_detect_complete(data[c("x", "ymin", "ymax")])
     ids <- cumsum(missing_pos) + 1
     ids[missing_pos] <- NA
 
diff --git a/R/ggunit.R b/R/ggunit.R
new file mode 100755
index 0000000000..30acd393e4
--- /dev/null
+++ b/R/ggunit.R
@@ -0,0 +1,296 @@
+# constructors ------------------------------------------------------------
+
+new_ggunit <- function(x = null_unit()) {
+  x <- vec_cast(x, null_unit())
+  class(x) <- c("ggunit", class(x), "vctrs_vctr")
+  x
+}
+
+ggunit <- function(x = numeric(), units = "native", data = NULL) {
+  x <- vec_cast(x, numeric())
+  units <- vec_cast(units, character())
+  data <- vec_cast(data, list())
+
+  if (length(x) == 0) {
+    x <- null_unit()
+  } else {
+    x <- unit(x, units, data = data)
+  }
+
+  new_ggunit(x)
+}
+
+
+# casting ----------------------------------------------------------------
+
+as_ggunit <- function(x) {
+  vec_cast(x, new_ggunit())
+}
+
+as_pt <- function(x) {
+  ggunit(x, "pt")
+}
+
+as_npc <- function(x) {
+  ggunit(x, "npc")
+}
+
+
+# type predicates ---------------------------------------------------------
+
+is_ggunit <- function(x) {
+  inherits(x, "ggunit")
+}
+
+
+# math --------------------------------------------------------------------
+
+#' @export
+Ops.ggunit <- function(x, y) {
+  if (!(.Generic %in% c("*", "/"))) {
+    x <- vec_cast(x, new_ggunit())
+    if (!missing(y)) {
+      y <- vec_cast(y, new_ggunit())
+    }
+  }
+
+  if (.Generic %in% c("==", ">", "<", ">=", "<=")) {
+    f <- match.fun(.Generic)
+    len <- max(length(x), length(y))
+    out <- logical(len)
+    x <- rep_len(x, len)
+    y <- rep_len(y, len)
+    type_x = unitType(x)
+    type_y = unitType(y)
+
+    is_same_type <- type_x == type_y & !type_x %in% c("sum", "min", "max")
+    out[is_same_type] <- f(as.numeric(x[is_same_type]), as.numeric(y[is_same_type]))
+
+    # determine relationships in otherwise incomparable units where possible (unequal signs, 0s, Infs)
+    not_same_type <- !is_same_type
+    x_not_same_type <- x[not_same_type]
+    y_not_same_type <- y[not_same_type]
+    sign_x_not_same_type <- sign(x_not_same_type)
+    sign_y_not_same_type <- sign(y_not_same_type)
+    out[not_same_type] <- switch(.Generic,
+      "==" =
+        sign_x_not_same_type == sign_y_not_same_type &
+        ifelse(sign_x_not_same_type == 0 | (is.infinite(x_not_same_type) & is.infinite(y_not_same_type)), TRUE, NA),
+      "<" =, ">" =
+        ifelse(sign_x_not_same_type == sign_y_not_same_type,
+          ifelse(sign_x_not_same_type == 0, FALSE, NA),
+          f(sign_x_not_same_type, sign_y_not_same_type)
+        ),
+      "<=" = x_not_same_type < y_not_same_type | x_not_same_type == y_not_same_type,
+      ">=" = x_not_same_type > y_not_same_type | x_not_same_type == y_not_same_type
+    )
+  } else {
+    out <- new_ggunit(NextMethod())
+  }
+
+  if (.Generic %in% c("+", "-")) {
+    out <- collapse_native_units(out)
+  }
+
+  out
+}
+
+#' @export
+chooseOpsMethod.ggunit = function(x, y, mx, my, cl, reverse) {
+  # TODO: something more comprehensive?
+  is_ggunit(vec_ptype2(x, y)) == is_ggunit(x)
+}
+
+#' @export
+Summary.ggunit <- function(..., na.rm = FALSE) {
+  ggunits <- vec_cast_common(..., .to = new_ggunit())
+  units <- vec_cast(ggunits, list_of(null_unit()))
+  args <- c(units, list(na.rm = na.rm))
+  out <- switch(.Generic,
+    range = vec_c(do.call(min, args), do.call(max, args)),
+    sum =, min =, max = collapse_native_units(do.call(.Generic, args)),
+    do.call(.Generic, args)
+  )
+  new_ggunit(out)
+}
+
+ggunit_math_function <- function(x, out_type, atomic_f, sum_f = function(x) NA, min_f = function(x) NA, max_f = function(x) NA) {
+  out <- out_type(length(x))
+  type <- unitType(x)
+
+  is_atomic <- !type %in% c("sum", "min", "max")
+  out[is_atomic] <- atomic_f(as.numeric(x[is_atomic]))
+
+  for (t in c("sum", "min", "max")) {
+    is_type <- t == type
+    f <- get(paste0(t, "_f"))
+    out[is_type] <- vapply(x[is_type], FUN.VALUE = out_type(1), function(x_i) {
+      components <- vec_cast(unclass(x_i)[[1]][[2]], new_ggunit())
+      f(components)
+    })
+  }
+
+  out
+}
+
+is_pos_Inf <- function(x) (is.infinite(x) & sign(x) == 1) %in% TRUE
+is_neg_Inf <- function(x) (is.infinite(x) & sign(x) == -1) %in% TRUE
+
+#' @export
+is.infinite.ggunit <- function(x) {
+  ggunit_math_function(x, logical,
+    atomic_f = is.infinite,
+    sum_f = function(x) any(is.infinite(x)),
+    min_f = function(x) all(is.infinite(x)) || any(is_neg_Inf(x)),
+    max_f = function(x) all(is.infinite(x)) || any(is_pos_Inf(x))
+  )
+}
+
+#' @export
+is.finite.ggunit <- function(x) {
+  ggunit_math_function(x, logical,
+    atomic_f = is.finite,
+    sum_f = function(x) all(is.finite(x)),
+    min_f = function(x) any(is.finite(x)) && all(is_pos_Inf(x[!is.finite(x)])),
+    max_f = function(x) any(is.finite(x)) && all(is_neg_Inf(x[!is.finite(x)]))
+  )
+}
+
+#' @export
+is.na.ggunit <- function(x) {
+  ggunit_math_function(x, logical, atomic_f = is.na, sum_f = anyNA, min_f = anyNA, max_f = anyNA)
+}
+
+#' @export
+sign.ggunit <- function(x) {
+  ggunit_math_function(x, numeric,
+    atomic_f = sign,
+    sum_f = function(x) {
+      unique_sign <- unique(sign(x))
+      if (length(unique_sign) == 1) unique_sign else NA_real_
+    },
+    min_f = function(x) {
+      sign_x <- sign(x)
+      if (isTRUE(any(sign_x == -1))) -1 else min(sign_x)
+    },
+    max_f = function(x) {
+      sign_x <- sign(x)
+      if (isTRUE(any(sign_x == 1))) 1 else max(sign_x)
+    }
+  )
+}
+
+ggunit_pmin <- function(...) {
+  dots <- vec_cast(list(...), list_of(new_ggunit()))
+  vec_cast(.mapply(min, dots, NULL), new_ggunit())
+}
+
+ggunit_pmax <- function(...) {
+  dots <- vec_cast(list(...), list_of(new_ggunit()))
+  vec_cast(.mapply(max, dots, NULL), new_ggunit())
+}
+
+
+# assignment --------------------------------------------------------------
+
+#' @export
+`[.ggunit` <- function(x, i) {
+  if (missing(i)) return(x)
+  vec_slice(x, i)
+}
+
+#' @export
+`[<-.ggunit` <- function(x, i, ..., value) {
+  value <- vec_cast(value, x)
+  out <- NextMethod()
+  new_ggunit(out)
+}
+
+#' @export
+`[[<-.ggunit` <- function(x, i, ..., value) {
+  value <- vec_cast(value, x)
+  out <- NextMethod()
+  new_ggunit(out)
+}
+
+
+# printing ----------------------------------------------------------------
+
+#' @export
+print.ggunit <- function(x, ...) {
+  # need to manually provide this rather than relying on print.vctrs_vctr()
+  # to bypass the printing method for grid::unit
+  obj_print(x, ...)
+  invisible(x)
+}
+
+
+# proxies -----------------------------------------------------------------
+
+#' @export
+vec_restore.ggunit <- function(x, ...) {
+  x <- NextMethod()
+  class(x) <- c("ggunit", class(x), "vctrs_vctr")
+  x
+}
+
+
+# casting -----------------------------------------------------------------
+
+#' @export
+vec_ptype2.ggunit.ggunit <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.ggunit.unit <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.unit.ggunit <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.ggunit.simpleUnit <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.simpleUnit.ggunit <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.ggunit.double <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.double.ggunit <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.ggunit.integer <- function(x, y, ...) new_ggunit()
+#' @export
+vec_ptype2.integer.ggunit <- function(x, y, ...) new_ggunit()
+
+#' @export
+vec_cast.ggunit.ggunit <- function(x, to, ...) x
+#' @export
+vec_cast.ggunit.unit <- function(x, to, ...) new_ggunit(x)
+#' @export
+vec_cast.unit.ggunit <- function(x, to, ...) `class<-`(x, setdiff(class(x), c("ggunit", "vctrs_vctr")))
+#' @export
+vec_cast.ggunit.simpleUnit <- function(x, to, ...) new_ggunit(x)
+#' @export
+vec_cast.simpleUnit.ggunit <- function(x, to, ...) `class<-`(x, setdiff(class(x), c("ggunit", "vctrs_vctr")))
+#' @export
+vec_cast.ggunit.integer <- function(x, to, ...) ggunit(x)
+#' @export
+vec_cast.integer.ggunit <- function(x, to, ...) as.integer(as.numeric(x))
+#' @export
+vec_cast.ggunit.double <- function(x, to, ...) ggunit(x)
+#' @export
+vec_cast.double.ggunit <- function(x, to, ...) as.numeric(x)
+#' @export
+vec_cast.ggunit.logical <- function(x, to, ...) ggunit(x)
+#' @export
+vec_cast.logical.ggunit <- function(x, to, ...) as.logical(as.numeric(x))
+#' @export
+vec_cast.ggunit.mapped_discrete <- function(x, to, ...) ggunit(x)
+#' @export
+vec_cast.mapped_discrete.ggunit <- function(x, to, ...) new_mapped_discrete(as.numeric(x))
+#' @export
+vec_cast.ggunit.character <- function(x, to, ...) stop_incompatible_cast(x, to, x_arg = "a", to_arg = "to")
+#' @export
+vec_cast.ggunit.list <- function(x, to, ...) {
+  is_na <- vapply(x, is.null, logical(1))
+  x[is_na] <- NA
+  x <- vec_cast(x, list_of(new_ggunit()))
+  if (any(lengths(x) != 1)) {
+    stop_incompatible_cast(x, to, x_arg = "x", to_arg = "to", details = "All elements of the list must be length-1 ggunits or NULL.")
+  }
+  list_unchop(x, ptype = new_ggunit())
+}
diff --git a/R/performance.R b/R/performance.R
index b26b1a7072..dff75d5de1 100644
--- a/R/performance.R
+++ b/R/performance.R
@@ -13,7 +13,7 @@ mat_2_df <- function(x, col_names = colnames(x)) {
 df_col <- function(x, name) .subset2(x, name)
 
 df_rows <- function(x, i) {
-  cols <- lapply(x, `[`, i = i)
+  cols <- lapply(x, vec_slice, i = i)
   data_frame0(!!!cols, .size = length(i))
 }
 
diff --git a/R/plot-build.R b/R/plot-build.R
index 700260c281..7b2f271f76 100644
--- a/R/plot-build.R
+++ b/R/plot-build.R
@@ -51,7 +51,7 @@ ggplot_build.ggplot <- function(plot) {
 
   # Compute aesthetics to produce data with generalised variable names
   data <- by_layer(function(l, d) l$compute_aesthetics(d, plot), layers, data, "computing aesthetics")
-  data <- .ignore_data(data)
+  data <- .ignore_data(.ignore_units(data))
 
   # Transform all scales
   data <- lapply(data, scales$transform_df)
@@ -63,7 +63,7 @@ ggplot_build.ggplot <- function(plot) {
 
   layout$train_position(data, scale_x(), scale_y())
   data <- layout$map_position(data)
-  data <- .expose_data(data)
+  data <- .expose_units(.expose_data(data))
 
   # Apply and map statistics
   data <- by_layer(function(l, d) l$compute_statistic(d, layout), layers, data, "computing stat")
@@ -76,6 +76,7 @@ ggplot_build.ggplot <- function(plot) {
   data <- by_layer(function(l, d) l$compute_geom_1(d), layers, data, "setting up geom")
 
   # Apply position adjustments
+  data <- .ignore_units(data)
   data <- by_layer(function(l, d) l$compute_position(d, layout), layers, data, "computing position")
 
   # Reset position scales, then re-train and map.  This ensures that facets
@@ -100,7 +101,7 @@ ggplot_build.ggplot <- function(plot) {
     # Only keep custom guides if there are no non-position scales
     plot$guides <- plot$guides$get_custom()
   }
-  data <- .expose_data(data)
+  data <- .expose_units(.expose_data(data))
 
   # Fill in defaults etc.
   data <- by_layer(function(l, d) l$compute_geom_2(d), layers, data, "setting up geom aesthetics")
diff --git a/R/position-.R b/R/position-.R
old mode 100644
new mode 100755
diff --git a/R/scale-.R b/R/scale-.R
old mode 100644
new mode 100755
diff --git a/R/scale-type.R b/R/scale-type.R
index 2feaa69c82..cffac38ec9 100644
--- a/R/scale-type.R
+++ b/R/scale-type.R
@@ -102,3 +102,6 @@ scale_type.double <- function(x) "continuous"
 
 #' @export
 scale_type.hms <- function(x) "time"
+
+#' @export
+scale_type.unit <- function(x) "continuous"
diff --git a/R/utilities-resolution.R b/R/utilities-resolution.R
index 28e54cd969..cf2486b83c 100644
--- a/R/utilities-resolution.R
+++ b/R/utilities-resolution.R
@@ -18,6 +18,9 @@
 #' resolution(c(2, 10, 20, 50))
 #' resolution(c(2L, 10L, 20L, 50L))
 resolution <- function(x, zero = TRUE) {
+  if (is.unit(x)) {
+    x <- native_units(x)
+  }
   if (is.integer(x) || is_mapped_discrete(x) ||
       zero_range(range(x, na.rm = TRUE))) {
     return(1)
diff --git a/R/utilities-unit.R b/R/utilities-unit.R
new file mode 100755
index 0000000000..c5d3e6041f
--- /dev/null
+++ b/R/utilities-unit.R
@@ -0,0 +1,216 @@
+#' transform x via the function trans. If x is a grid::unit(), apply the
+#' transformation only to "native" units within x.
+#' @noRd
+transform_native_units <- function(x, trans, ...) {
+  if (!is.unit(x)) {
+    return(trans(x, ...))
+  }
+
+  x <- collapse_native_units(x)
+  native_units(x) <- trans(native_units(x), ...)
+  x
+}
+
+unit_components <- function(x) {
+  unclass(x)[[1]][[2]]
+}
+
+`unit_components<-` <- function(x, value) {
+  # force the value to be a list form of unit, not a simpleUnit
+  x <- vec_cast(x, null_unit())
+
+  oldclass <- class(x)
+  x <- unclass(x)
+  x[[1]][[2]] <- value
+  class(x) <- oldclass
+  x
+}
+
+collapse_native_units <- function(x) {
+  x <- vec_cast(x, new_ggunit())
+  type <- unitType(x)
+  is_recursive <- type %in% c("sum", "min", "max")
+  x[is_recursive] <- .mapply(list(x[is_recursive], type[is_recursive]), NULL, FUN = function(x_i, f) {
+    f <- match.fun(f)
+    components <- collapse_native_units(unit_components(x_i))
+    is_native <- unitType(components) == "native"
+    if (any(is_native)) {
+      unit_components(x_i) <- vec_c(unit(f(as.numeric(components[is_native])), "native"), components[!is_native])
+    }
+    x_i
+  })
+  x
+}
+
+native_units <- function(x) {
+  if (!is.unit(x)) {
+    if (is.numeric(x)) return(x)
+    stop_input_type(x, as_cli("a {.cls unit} or a {.cls numeric}"))
+  }
+  .get_native_units(x)$values
+}
+
+.get_native_units <- function(x) {
+  values <- rep_len(NA_real_, length(x))
+  type <- unitType(x)
+
+  is_native <- type == "native"
+  if (any(is_native)) {
+    values[is_native] <- as.numeric(x[is_native])
+  }
+
+  is_recursive <- unitType(x) %in% c("sum", "min", "max")
+  if (any(is_recursive)) {
+    for (i in which(is_recursive)) {
+      out <- .get_native_units(unit_components(x[[i]]))
+      native_i <- which(out$is_native)
+      if (length(native_i) > 1) {
+        cli::cli_abort("More than one native unit in {x[[i]]}")
+      } else if (length(native_i) == 1) {
+        values[[i]] <- out$values[[native_i]]
+        is_native[[i]] <- TRUE
+      }
+    }
+  }
+
+  list(values = values, is_native = is_native)
+}
+
+`native_units<-` <- function(x, values) {
+  if (!is.unit(x)) {
+    if (is.numeric(x)) return(values)
+    stop_input_type(x, as_cli("a {.cls unit} or a {.cls numeric}"))
+  }
+  .set_native_units(x, values)$x
+}
+
+.set_native_units <- function(x, values) {
+  len <- max(length(x), length(values))
+  x <- rep_len(x, len)
+  values <- rep_len(values, len)
+  type <- unitType(x)
+
+  is_native <- type == "native"
+  if (any(is_native)) {
+    x[is_native] <- unit(values[is_native], "native")
+  }
+
+  is_recursive <- unitType(x) %in% c("sum", "min", "max")
+  if (any(is_recursive)) {
+    for (i in which(is_recursive)) {
+      out <- .set_native_units(unit_components(x[[i]]), values[[i]])
+      native_i <- which(out$is_native)
+      if (length(native_i) > 1) {
+        cli::cli_abort("More than one native unit in {x[[i]]}")
+      } else if (length(native_i) == 1) {
+        unit_components(x[[i]]) <- out$x
+        is_native[[i]] <- TRUE
+      }
+    }
+  }
+
+  list(x = x, is_native = is_native)
+}
+
+.ignore_units <- function(data, cols = c(ggplot_global$x_aes, ggplot_global$y_aes)) {
+  if (is.data.frame(data)) {
+    return(.ignore_units(list(data), cols)[[1]])
+  }
+
+  lapply(data, function(df) {
+    if (is.null(cols)) {
+      is_selected <- TRUE
+    } else {
+      is_selected <- names(df) %in% cols
+    }
+    is_unit <- vapply(df, is.unit, logical(1)) & is_selected
+    if (!any(is_unit)) {
+      return(df)
+    }
+    df <- unclass(df)
+    unit_cols <- lapply(df[is_unit], collapse_native_units)
+    new_data_frame(c(
+      df[!is_unit],
+      lapply(unit_cols, native_units),
+      list(.ignored_units = new_data_frame(unit_cols))
+    ))
+  })
+}
+
+.expose_units <- function(data) {
+  if (is.data.frame(data)) {
+    return(.expose_units(list(data))[[1]])
+  }
+
+  lapply(data, function(df) {
+    is_ignored <- which(names(df) == ".ignored_units")
+    if (length(is_ignored) == 0) {
+      return(df)
+    }
+    unit_col_names <- intersect(names(df), names(df[[is_ignored[1]]]))
+    is_unit <- which(names(df) %in% unit_col_names)
+    df <- unclass(df)
+    new_data_frame(c(
+      df[-c(is_ignored, is_unit)],
+      mapply(`native_units<-`, df[[is_ignored[1]]][unit_col_names], df[unit_col_names], SIMPLIFY = FALSE)
+    ))
+  })
+}
+
+
+# rescale -----------------------------------------------------------------
+
+#' @export
+rescale.unit <- function(x, to, from, ...) {
+  native_units(x) <- rescale(native_units(x), to, from, ...)
+  x
+}
+
+# proxies -----------------------------------------------------------------
+
+#' @export
+vec_proxy.unit <- function(x, ...) {
+  unclass(x)
+}
+
+#' @export
+vec_restore.unit <- function(x, ...) {
+  # replace NAs (NULL entries) with unit's version of NA
+  is_na <- vapply(x, is.null, logical(1))
+  x[is_na] <- vec_proxy(unit(NA_real_, "native"))
+
+  class(x) <- c("unit", "unit_v2")
+  x
+}
+
+#' @export
+vec_proxy.simpleUnit <- function(x, ...) {
+  # turn a simpleUnit into a unit when proxied, because simpleUnit's format
+  # (a numeric vector with an attribute indicating the type of all entries)
+  # does not work properly with many operations, like binding
+  type <- attr(x, "unit")
+  lapply(unclass(x), function(x_i) list(x_i, NULL, type))
+}
+
+
+# casting -----------------------------------------------------------------
+
+null_unit <- function() {
+  # grid::unit() doesn't allow zero-length vectors,
+  # so we have to do this manually
+  structure(list(), class = c("unit", "unit_v2"))
+}
+
+#' @export
+vec_ptype2.unit.unit <- function(x, y, ...) null_unit()
+#' @export
+vec_ptype2.unit.simpleUnit <- function(x, y, ...) null_unit()
+#' @export
+vec_ptype2.simpleUnit.unit <- function(x, y, ...) null_unit()
+
+#' @export
+vec_cast.unit.unit <- function(x, to, ...) x
+#' @export
+vec_cast.unit.simpleUnit <- function(x, to, ...) vec_restore(vec_proxy(x), null_unit())
+#' @export
+vec_cast.simpleUnit.unit <- function(x, to, ...) vec_restore(vec_proxy(x), null_unit())