diff --git a/NAMESPACE b/NAMESPACE
index 3c2cf0128..87e852406 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,6 +9,7 @@ S3method("[",sf)
 S3method("[",sfc)
 S3method("[<-",sf)
 S3method("[<-",sfc)
+S3method("[[",sfc)
 S3method("[[<-",sf)
 S3method("st_agr<-",sf)
 S3method("st_crs<-",bbox)
diff --git a/R/arith.R b/R/arith.R
index e6df6ed5a..721eed037 100644
--- a/R/arith.R
+++ b/R/arith.R
@@ -133,6 +133,11 @@ Ops.sfc <- function(e1, e2) {
 	if (length(e1) == 0) # empty set
 		return(e1)
 
+	if (inherits(e1, "sfc")) # realize
+		e1 = e1[]
+	if (inherits(e2, "sfc"))
+		e2 = e2[]
+
 	if (is.numeric(e2) && !is.matrix(e2) && length(e2) <= 2 && .Generic %in% c("+", "-")) {
 		if (.Generic == "-")
 			e2 <- -e2
diff --git a/R/bbox.R b/R/bbox.R
index 69e265ea6..a91697398 100644
--- a/R/bbox.R
+++ b/R/bbox.R
@@ -9,6 +9,13 @@ bb_wrap = function(bb) {
 	structure(as.double(bb), names = c("xmin", "ymin", "xmax", "ymax"), class = "bbox")
 }
 
+bbox.pointmatrix = function(obj, ...) {
+	if (nrow(obj) == 0 || all(is.na(obj)))
+		bb_wrap(rep(NA_real_, 4))
+	else
+		bb_wrap(as.vector(t(apply(obj[,1:2,drop=FALSE], 2, range, na.rm = TRUE))))
+}
+
 bbox.Set = function(obj, ...) {
 	sel = !sfc_is_empty(obj)
 	if (! any(sel))
@@ -134,7 +141,9 @@ print.bbox = function(x, ...) {
 }
 
 compute_bbox = function(obj) {
-	switch(class(obj)[1],
+	if (!is.null(pts <- attr(obj, "points")))
+		bbox.pointmatrix(pts)
+	else switch(class(obj)[1],
 		sfc_POINT = bb_wrap(bbox.Set(obj)),
 		sfc_MULTIPOINT = bb_wrap(bbox.MtrxSet(obj)),
 		sfc_LINESTRING = bb_wrap(bbox.MtrxSet(obj)),
diff --git a/R/bind.R b/R/bind.R
index f5fc22448..a7fab38f4 100644
--- a/R/bind.R
+++ b/R/bind.R
@@ -36,7 +36,12 @@ rbind.sf = function(..., deparse.level = 1) {
 		else
 			NULL
 	chk_equal_crs(dots)
-	ret = st_sf(rbind.data.frame(...), crs = st_crs(dots[[1]]), sf_column_name = sf_column)
+	crs0 = st_crs(dots[[1]])
+	for (i in seq_along(dots)) {
+		if (all(sapply(unclass(st_geometry(dots[[i]])), is.null)))
+			st_geometry(dots[[i]]) = st_geometry(dots[[i]])[] # realize
+	}
+	ret = st_sf(do.call(rbind.data.frame, dots), crs = crs0, sf_column_name = sf_column)
 	st_geometry(ret) = st_sfc(st_geometry(ret)) # might need to reclass to GEOMETRY
 	bb = do.call(rbind, lapply(dots, st_bbox))
 	bb = bb_wrap(c(min(bb[,1L], na.rm = TRUE), min(bb[,2L], na.rm = TRUE),
diff --git a/R/geom-measures.R b/R/geom-measures.R
index 4c85df450..b353215f5 100644
--- a/R/geom-measures.R
+++ b/R/geom-measures.R
@@ -206,6 +206,7 @@ st_distance = function(x, y, ..., dist_fun, by_element = FALSE,
 				else
 					CPL_geos_dist(x, y, which, par)
 			}
+		d[is.nan(d)] = NA_real_
 		if (!is.null(u <- st_crs(x)$ud_unit))
 			units(d) = u
 		d
@@ -254,7 +255,7 @@ st_line_project = function(line, point, normalized = FALSE) {
 		is.logical(normalized), length(normalized) == 1,
 		st_crs(line) == st_crs(point))
 	line = st_cast(line, "LINESTRING")
-	point = st_cast(point, "POINT")
+	point = st_cast(point[], "POINT")
 	if (isTRUE(st_is_longlat(line)))
 		message_longlat("st_project_point")
 	recycled = recycle_common(list(line, point))
diff --git a/R/geom-transformers.R b/R/geom-transformers.R
index adb238bb7..24475ecb0 100644
--- a/R/geom-transformers.R
+++ b/R/geom-transformers.R
@@ -473,7 +473,7 @@ st_minimum_rotated_rectangle.sf = function(x, dTolerance, ...) {
 #'  n = 100
 #'  pts = st_as_sf(data.frame(matrix(runif(n), , 2), id = 1:(n/2)), coords = c("X1", "X2"))
 #'  # compute Voronoi polygons:
-#'  pols = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts))))
+#'  pols = st_collection_extract(st_voronoi(st_combine(pts)))
 #'  # match them to points:
 #'  pts$pols = pols[unlist(st_intersects(pts, pols))]
 #'  plot(pts["id"], pch = 16) # ID is color
@@ -744,8 +744,14 @@ st_segmentize.sf = function(x, dfMaxLength, ...) {
 #' @examples
 #' nc = st_read(system.file("shape/nc.shp", package="sf"))
 #' st_combine(nc)
-st_combine = function(x)
-	st_sfc(do.call(c, st_geometry(x)), crs = st_crs(x)) # flatten/merge
+st_combine = function(x) {
+	x = st_geometry(x)
+	if (inherits(x, "sfc_POINT") && !is.null(pts <- attr(x, "points"))) {
+		mp = structure(list(st_multipoint(pts, attr(x, "point_dim"))), bbox = st_bbox(x))
+		st_sfc(mp, crs = st_crs(x))
+	} else
+		st_sfc(do.call(c, x), crs = st_crs(x)) # flatten/merge
+}
 
 # x: object of class sf
 # y: object of class sf or sfc
diff --git a/R/m_range.R b/R/m_range.R
index e7c71bba6..4912a89da 100644
--- a/R/m_range.R
+++ b/R/m_range.R
@@ -9,6 +9,10 @@ mb_wrap = function(mb) {
 	structure(mb, names = c("mmin", "mmax"), class = "m_range")
 }
 
+m_range.pointmatrix = function(obj, ...) {
+	mb_wrap(range(obj[,ncol(obj)]))
+}
+
 m_range.Set = function(obj, ...) {
 	sel = vapply(obj, function(x) { length(x) && !all(is.na(x)) }, TRUE)
 	if (! any(sel))
@@ -130,7 +134,9 @@ print.m_range = function(x, ...) {
 }
 
 compute_m_range = function(obj) {
-	switch(class(obj)[1],
+	if (!is.null(pts <- attr(obj, "points")))
+		m_range.pointmatrix(pts)
+	else switch(class(obj)[1],
 		   sfc_POINT = mb_wrap(m_range.Set(obj)),
 		   sfc_MULTIPOINT = mb_wrap(m_range.MtrxSet(obj)),
 		   sfc_LINESTRING = mb_wrap(m_range.MtrxSet(obj)),
diff --git a/R/plot.R b/R/plot.R
index 32cecc78e..aeec2ea14 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -345,7 +345,8 @@ plot.sfc_POINT = function(x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1,
 	col = rep(col, length.out = npts)
 	bg = rep(bg, length.out = npts)
 	cex = rep(cex, length.out = npts)
-	mat = t(matrix(unlist(x, use.names = FALSE), ncol = length(x))) #933
+	#mat = t(matrix(unlist(x, use.names = FALSE), ncol = length(x))) #933
+	mat = st_coordinates(x)
 	if (!is.null(mat)) {
 		ne = !is.na(rowMeans(mat))  ## faster than apply; #933
 		points(mat[ne,, drop = FALSE], pch = pch[ne], col = col[ne], bg = bg[ne],
diff --git a/R/sf.R b/R/sf.R
index eb262dced..96db45d27 100644
--- a/R/sf.R
+++ b/R/sf.R
@@ -46,24 +46,31 @@ st_as_sf.data.frame = function(x, ..., agr = NA_agr_, coords, wkt,
 		else
 			x$geometry = st_as_sfc(as.character(x[[wkt]]))
 	} else if (! missing(coords)) {
-		cc = as.data.frame(lapply(x[coords], as.numeric))
-		if (na.fail && anyNA(cc))
-			stop("missing values in coordinates not allowed")
-		# classdim = getClassDim(rep(0, length(coords)), length(coords), dim, "POINT")
+		if (length(coords) == 1) {
+			stopifnot(is.matrix(x[[coords]]), is.numeric(x[[coords]]))
+			cc = x[[coords]]
+		} else {
+			if (length(coords) == 2)
+				dim = "XY"
+			stopifnot(length(coords) == nchar(dim), dim %in% c("XY", "XYZ", "XYZM", "XYM"))
+			cc = do.call(cbind, lapply(x[coords], as.numeric))
+			if (na.fail && anyNA(cc))
+				stop("missing values in coordinates not allowed")
+		}
+		dimnames(cc) = NULL
 		if (is.null(sf_column_name))
 			sf_column_name = "geometry"
+
 		x[[sf_column_name]] = if (nchar(dim) < 4 && ncol(cc) == 4) { # create POLYGONs:
 				fn = function(x) st_as_sfc(st_bbox(c(xmin = x[[1]], ymin = x[[2]], xmax = x[[3]], ymax = x[[4]])))
 				do.call(c, apply(as.matrix(cc), 1, fn))
 			} else { # points:
-				structure( points_rcpp(as.matrix(cc), dim),
+				structure(vector("list", length = nrow(cc)),
+				points = cc,
+				points_dim = dim,
 				n_empty = 0L, precision = 0, crs = NA_crs_,
-				bbox = structure(
-					c(xmin = min(cc[[1]], na.rm = TRUE),
-					ymin = min(cc[[2]], na.rm = TRUE),
-					xmax = max(cc[[1]], na.rm = TRUE),
-					ymax = max(cc[[2]], na.rm = TRUE)), class = "bbox"),
-				class =  c("sfc_POINT", "sfc" ), names = NULL)
+				bbox = bbox.pointmatrix(cc),
+				class = c("sfc_POINT", "sfc"), names = NULL)
 			}
 
 		if (remove) {
@@ -430,9 +437,10 @@ print.sf = function(x, ..., n = getOption("sf_max_print", default = 10)) {
 		app = paste0(app, "\n", "Active geometry column: ", attr(x, "sf_column"))
 	print(st_geometry(x), n = 0, what = "Simple feature collection with", append = app)
 	if (n > 0) {
-		if (inherits(x, c("tbl_df", "tbl")))
+		if (inherits(x, c("tbl_df", "tbl"))) {
+			st_geometry(x) = x[[attr(x, "sf_column")]][] # note the extra []: this reloads points
 			NextMethod()
-		else {
+		} else {
 			y <- x
 			if (nrow(y) > n) {
 				cat(paste("First", n, "features:\n"))
diff --git a/R/sfc.R b/R/sfc.R
index c3230fb2a..7f6346cde 100644
--- a/R/sfc.R
+++ b/R/sfc.R
@@ -54,14 +54,24 @@ st_sfc = function(..., crs = NA_crs_, precision = 0.0, check_ring_dir = FALSE, d
 		lst = lst[[1]]
 	stopifnot(is.numeric(crs) || is.character(crs) || inherits(crs, "crs"))
 
+	points_in_attr <- !is.null(attr(lst, "points"))
+
 	# check for NULLs:
 	a = attributes(lst)
-	is_null = sfc_is_null(lst)
+	is_null = if (points_in_attr)
+			rep(FALSE, length(lst))
+		else
+			sfc_is_null(lst)
 	lst = unclass(lst)
-	lst = lst[! is_null]
+	if (!points_in_attr)
+		lst = lst[! is_null]
+
 	attributes(lst) = a
 
-	dims_and_types = sfc_unique_sfg_dims_and_types(lst)
+	dims_and_types = if (points_in_attr)
+			list(class_dim = attr(lst, "points_dim"), class_type = "POINT")
+		else
+			sfc_unique_sfg_dims_and_types(lst)
 	
 	cls = if (length(lst) == 0) # empty set: no geometries read
 		c("sfc_GEOMETRY", "sfc")
@@ -155,21 +165,46 @@ st_sfc = function(..., crs = NA_crs_, precision = 0.0, check_ring_dir = FALSE, d
 #' @details if `x` has a `dim` attribute (i.e. is an `array` or `matrix`) then `op` cannot be used.
 #' @export
 "[.sfc" = function(x, i, j, ..., op = st_intersects) {
-	precision = st_precision(x)
-	crs = st_crs(x)
-	dim = if (length(x)) class(x[[1]])[1] else "XY"
-	if (!missing(i) && (inherits(i, c("sf", "sfc", "sfg"))))
+
+	if (!missing(i) && (inherits(i, "sf") || inherits(i, "sfc") || inherits(i, "sfg")))
 		i = lengths(op(x, i, ...)) != 0
-	if (!is.null(dim(x))) # x is an array with geometries
-		st_sfc(NextMethod(), crs = crs, precision = precision, dim = dim)
-	else # x is a list but avoid NextMethod() to allow j, ... to be specified & ignored:
-		st_sfc(unclass(x)[i], crs = crs, precision = precision, dim = dim)
+	if (inherits(x, "sfc_POINT") && !is.null(attr(x, "points")))
+		st_sfc(restore_points(x, i), crs = st_crs(x), precision = st_precision(x),
+			dim = if(length(x)) class(x[[1]])[1] else "XY")
+	else {
+		precision = st_precision(x)
+		crs = st_crs(x)
+		dim = if (length(x)) class(x[[1]])[1] else "XY"
+		if (!is.null(dim(x))) # x is an array with geometries
+			st_sfc(NextMethod(), crs = crs, precision = precision, dim = dim)
+		else # x is a list but avoid NextMethod() to allow j, ... to be specified & ignored:
+			st_sfc(unclass(x)[i], crs = crs, precision = precision, dim = dim)
+	}
 }
 
 #' @export
-"[<-.sfc" = function (x, i, j, ..., value) {
-	if (is.null(value) || inherits(value, "sfg"))
+"[<-.sfc" = function (x, i, value) {
+	if (is.null(value) || inherits(value, "sfg")) {
+		is_points = inherits(value, "POINT")
 		value = list(value)
+	} else
+		is_points = inherits(value, "sfc_POINT")
+	if (inherits(x, "sfc_POINT") && !is.null(attr(x, "points"))) {
+		if (is_points) {
+			repl = if (!is.null(pts <- attr(value, "points")))
+					pts
+				else
+					do.call(rbind, value)
+			attr(x, "points")[i, ] = repl
+			return(structure(x, 
+							 n_empty = sum(is.na(attr(x, "points")[,1])), 
+							 bbox = bbox.pointmatrix(attr(x, "points"))
+							 )) # RETURNS
+		} else
+			x = x[] # realize
+	} 
+	value = value[] # realize in case sfc_POINT while x is not 
+	x = unclass(x) # becomes a list, but keeps attributes
 	ret = st_sfc(NextMethod(), recompute_bbox = TRUE)
 	structure(ret, n_empty = sum(sfc_is_empty(ret)))
 }
@@ -188,8 +223,17 @@ c.sfc = function(..., recursive = FALSE) {
 		else
 			c(ucls, "sfc")
 
+	points_attr = sapply(lst, function(x) !is.null(attr(x, "points")))
+	if (any(points_attr) && !all(points_attr)) {
+		for (i in seq_along(lst))
+			lst[[i]] = lst[[i]][] # realize
+		points_attr = FALSE
+	}
+
 	ret = unlist(lapply(lst, unclass), recursive = FALSE)
 	attributes(ret) = attributes(lst[[1]]) # crs
+	if (all(points_attr))
+		attr(ret, "points") = do.call(rbind, lapply(lst, attr, "points"))
 	class(ret) = cls
 	attr(ret, "bbox") = compute_bbox(ret) # dispatch on class
 	attr(ret, "n_empty") = sum(sapply(lst, attr, which = "n_empty"))
@@ -518,8 +562,11 @@ st_coordinates.sfc = function(x, ...) {
 		return(matrix(nrow = 0, ncol = 2))
 
 	ret = switch(class(x)[1],
-		sfc_POINT = matrix(unlist(x, use.names = FALSE), nrow = length(x), byrow = TRUE,
-		     dimnames = NULL),
+		sfc_POINT = if (is.null(attr(x, "points"))) {
+					matrix(unlist(x, use.names = FALSE), nrow = length(x), byrow = TRUE, dimnames = NULL)
+				} else {
+					attr(x, "points")
+				},
 		sfc_MULTIPOINT = ,
 		sfc_LINESTRING = coord_2(x),
 		sfc_MULTILINESTRING = ,
@@ -527,7 +574,10 @@ st_coordinates.sfc = function(x, ...) {
 		sfc_MULTIPOLYGON = coord_4(x),
 		stop(paste("not implemented for objects of class", class(x)[1]))
 	)
-	Dims = class(x[[1]])[1]
+	Dims = if (!is.null(attr(x, "points_dim")))
+			attr(x, "points_dim")
+		else
+			class(x[[1]])[1]
 	ncd = nchar(Dims)
 	colnames(ret)[1:ncd] = vapply(seq_len(ncd), function(i) substr(Dims, i, i), "")
 	ret
@@ -663,3 +713,24 @@ st_is_full.sf = function(x, ...) {
 st_is_full.bbox = function(x, ...) {
 	sf_use_s2() && st_is_longlat(x) && all(x == c(-180,-90,180,90))
 }
+
+#' @export
+`[[.sfc` = function(x, i, j, ..., exaxt = TRUE) {
+	if (inherits(x, "sfc_POINT") && !is.null(attr(x, "points")))
+		restore_point(x, i)
+	else
+		NextMethod()
+}
+
+restore_point = function(x, i = TRUE) {
+	restore_points(x, i)[[1]]
+}
+
+restore_points = function(x, i = TRUE) {
+	a = attributes(x)
+	points = a$points[i, , drop=FALSE]
+	structure(points_rcpp(points, a$points_dim),
+		n_empty = 0L, precision = a$precision, crs = a$crs,
+		bbox = bbox.pointmatrix(points), class = a$class, 
+		points = NULL, points_dim = NULL)
+}
diff --git a/R/sfg.R b/R/sfg.R
index de77e3f72..666995d08 100644
--- a/R/sfg.R
+++ b/R/sfg.R
@@ -258,6 +258,8 @@ c.sfg = function(..., recursive = FALSE, flatten = TRUE) {
 	Paste0 = function(lst) lapply(lst, unclass)
 	Paste1 = function(lst) do.call(c, lapply(lst, unclass))
 	lst = list(...)
+	if (length(lst) && is.null(lst[[1]]))
+		stop("to combine POINTs into MULTIPOINT, use st_combine(), or realize them first using x[]")
 	if (flatten) {
 		cls = vapply(lst, function(x) class(x)[2], "")
 		ucls = unique(cls)
diff --git a/R/sp.R b/R/sp.R
index 03ed81a5f..cd37ebca1 100644
--- a/R/sp.R
+++ b/R/sp.R
@@ -93,9 +93,10 @@ handle_bbox = function(sfc, sp) {
 st_as_sfc.SpatialPoints = function(x, ..., precision = 0.0) {
 	cc = x@coords
 	dimnames(cc) = NULL
-	lst = lapply(seq_len(nrow(cc)), function(x) st_point(cc[x,]))
-	handle_bbox(do.call(st_sfc, append(lst, list(crs = st_crs(x@proj4string), 
-		precision = precision))), x)
+	lst = vector("list", length = nrow(cc))
+	attr(lst, "points") = cc
+	attr(lst, "points_dim") = "XY"
+	handle_bbox(structure(st_sfc(lst, crs = st_crs(x@proj4string), precision = precision), n_empty = 0), x)
 }
 
 #' @rdname st_as_sfc
@@ -341,7 +342,11 @@ as_Spatial = function(from, cast = TRUE, IDs = paste0("ID", seq_along(from))) {
 sfc2SpatialPoints = function(from, IDs) {
 	if (!requireNamespace("sp", quietly = TRUE))
 		stop("package sp required, please install it first")
-	sp::SpatialPoints(do.call(rbind, from), proj4string = as(st_crs(from), "CRS"))
+	m = if (!is.null(pts <- attr(from, "points")))
+			pts
+		else 
+			do.call(rbind, from)
+	sp::SpatialPoints(m, proj4string = as(st_crs(from), "CRS"))
 }
 
 sfc2SpatialMultiPoints = function(from) {
diff --git a/R/tidyverse-vctrs.R b/R/tidyverse-vctrs.R
index 23509ec65..3d7b6fb85 100644
--- a/R/tidyverse-vctrs.R
+++ b/R/tidyverse-vctrs.R
@@ -3,14 +3,13 @@
 # time, this declares `sfc` lists as vectors which is necessary
 # because vctrs generally treats S3 lists as scalars.
 vec_proxy.sfc = function(x, ...) {
-	x
+	x[]
 }
 # This restores `sfc` attributes after manipulation of the proxy
 # (e.g. slicing or combination)
 vec_restore.sfc = function(x, to, ...) {
 	# Ensure restoration of `n_empty` by `st_sfc()`
 	attr(x, "n_empty") = NULL
-
 	st_sfc(x, crs = st_crs(to), precision = st_precision(to))
 }
 
diff --git a/R/tidyverse.R b/R/tidyverse.R
index 7e2ddc83d..9cfdf83f7 100644
--- a/R/tidyverse.R
+++ b/R/tidyverse.R
@@ -43,8 +43,14 @@ dplyr_reconstruct.sf = function(data, template) {
 #' }
 filter.sf <- function(.data, ..., .dots) {
 	agr = st_agr(.data)
+	g = st_geometry(.data)
 	class(.data) <- setdiff(class(.data), "sf")
-	.re_sf(NextMethod(), sf_column_name = attr(.data, "sf_column"), agr)
+	if (inherits(g, "sfc_POINT") && !is.null(pts <- attr(g, "points"))) {
+		.data[[ attr(.data, "sf_column") ]] = pts
+		st_as_sf(NextMethod(), coords = attr(.data, "sf_column"), agr = agr, remove = FALSE,
+			crs = st_crs(g)) # FIXME: doesn't handle tibble?
+	} else
+		.re_sf(NextMethod(), sf_column_name = attr(.data, "sf_column"), agr)
 }
 
 #' @name tidyverse
@@ -561,6 +567,8 @@ nest.sf = function(.data, ...) {
 	if (!requireNamespace("tidyr", quietly = TRUE))
 		stop("tidyr required: install first?")
 
+	if (inherits(g <- st_geometry(.data), "sfc_POINT") && !is.null(attr(g, "points")))
+		st_geometry(.data) = g[] # realize
 	class(.data) <- setdiff(class(.data), "sf")
 	ret = tidyr::nest(.data, ...)
 	lst = which(sapply(ret, inherits, "list"))[1]
@@ -650,7 +658,7 @@ type_sum.sfc <- function(x, ...) {
 
 #' @rdname tibble
 obj_sum.sfc <- function(x) {
-	vapply(x, function(sfg) format(sfg, width = 15L), "")
+	vapply(x[], function(sfg) format(sfg, width = 15L), "")
 }
 
 #' @rdname tibble
diff --git a/R/z_range.R b/R/z_range.R
index 23d8b1ffa..4724370f0 100644
--- a/R/z_range.R
+++ b/R/z_range.R
@@ -9,6 +9,10 @@ zb_wrap = function(zb) {
 	structure(zb, names = c("zmin", "zmax"), class = "z_range")
 }
 
+z_range.pointmatrix = function(obj, ...) {
+	zb_wrap(range(obj[,3]))
+}
+
 z_range.Set = function(obj, ...) {
 	sel = vapply(obj, function(x) { length(x) && !all(is.na(x)) }, TRUE)
 	if (! any(sel))
@@ -129,7 +133,9 @@ print.z_range = function(x, ...) {
 }
 
 compute_z_range = function(obj) {
-	switch(class(obj)[1],
+	if (!is.null(pts <- attr(obj, "points")))
+		z_range.pointmatrix(pts)
+	else switch(class(obj)[1],
 		   sfc_POINT = zb_wrap(z_range.Set(obj)),
 		   sfc_MULTIPOINT = zb_wrap(z_range.MtrxSet(obj)),
 		   sfc_LINESTRING = zb_wrap(z_range.MtrxSet(obj)),
diff --git a/man/geos_unary.Rd b/man/geos_unary.Rd
index 0041b1510..b1bccc425 100644
--- a/man/geos_unary.Rd
+++ b/man/geos_unary.Rd
@@ -276,7 +276,7 @@ if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.5.0") > -1) {
  n = 100
  pts = st_as_sf(data.frame(matrix(runif(n), , 2), id = 1:(n/2)), coords = c("X1", "X2"))
  # compute Voronoi polygons:
- pols = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts))))
+ pols = st_collection_extract(st_voronoi(st_combine(pts)))
  # match them to points:
  pts$pols = pols[unlist(st_intersects(pts, pols))]
  plot(pts["id"], pch = 16) # ID is color
diff --git a/src/wkb.cpp b/src/wkb.cpp
index c91277ff2..09debef7a 100644
--- a/src/wkb.cpp
+++ b/src/wkb.cpp
@@ -430,7 +430,32 @@ Rcpp::List CPL_read_wkb(Rcpp::List wkb_list, bool EWKB = false, bool spatialite
 	output.attr("n_empty") = (int) n_empty;
 	if ((EWKB || spatialite) && srid != 0)
 		output.attr("srid") = (int) srid;
-	return output;
+	if (n_types <= 1 && type == SF_Point) { // xxx
+		using namespace Rcpp; // so that later on the (i,_) works
+		NumericVector pt = output[0];
+		CharacterVector cls = pt.attr("class");
+		int nc = 2;
+		if (cls[0] == "XYZ" || cls[0] == "XYM")
+			nc = 3;
+		if (cls[0] == "XYZM")
+			nc = 4;
+		NumericMatrix m(wkb_list.size(), nc);
+		for (int i = 0; i < wkb_list.size(); i++) {
+			NumericVector nv = output[i];
+			m(i, _) = nv;
+		}
+		Rcpp::List ret(wkb_list.size());
+		ret.attr("points") = m;
+		CharacterVector points_dim(1);
+		points_dim(0) = cls(0);
+		ret.attr("points_dim") = points_dim;
+		ret.attr("single_type") = n_types <= 1; // if 0, we have only empty geometrycollections
+		ret.attr("n_empty") = (int) n_empty;
+		if ((EWKB || spatialite) && srid != 0)
+			ret.attr("srid") = (int) srid;
+		return ret; 
+	} else 
+		return output;
 }
 
 //
@@ -700,18 +725,39 @@ Rcpp::List CPL_write_wkb(Rcpp::List sfc, bool EWKB = false) {
 		}
 	}
 
-	for (int i = 0; i < sfc.size(); i++) {
-		Rcpp::checkUserInterrupt();
-		std::ostringstream os;
-		if (have_classes)
-			cls = classes[i];
-		write_data(os, sfc, i, EWKB, endian, cls, dm, precision, srid);
-		Rcpp::RawVector raw(os.str().size()); // os -> raw:
-		std::string str = os.str();
-		const char *cp = str.c_str();
-		for (size_t j = 0; j < str.size(); j++)
-			raw[j] = cp[j];
-		output[i] = raw; // raw vector to list
+	if (sfc.attr("points") != R_NilValue) {
+		using namespace Rcpp; // so that later on the (i,_) works
+		Rcpp::NumericMatrix m = sfc.attr("points");
+		unsigned int sf_type = make_type(cls, dm, EWKB, NULL, srid);
+		for (int i = 0; i < m.nrow(); i++) {
+			std::ostringstream os;
+			add_byte(os, (char) endian);
+			add_int(os, sf_type);
+			if (EWKB && srid != 0)
+				add_int(os, srid);
+			for (int j = 0; j < m.ncol(); j++)
+				add_double(os, m(i, j), precision);
+			Rcpp::RawVector raw(os.str().size()); // os -> raw:
+			std::string str = os.str();
+			const char *cp = str.c_str();
+			for (size_t j = 0; j < str.size(); j++)
+				raw[j] = cp[j];
+			output[i] = raw; // raw vector to list
+		}
+	} else {
+		for (int i = 0; i < sfc.size(); i++) {
+			Rcpp::checkUserInterrupt();
+			std::ostringstream os;
+			if (have_classes)
+				cls = classes[i];
+			write_data(os, sfc, i, EWKB, endian, cls, dm, precision, srid);
+			Rcpp::RawVector raw(os.str().size()); // os -> raw:
+			std::string str = os.str();
+			const char *cp = str.c_str();
+			for (size_t j = 0; j < str.size(); j++)
+				raw[j] = cp[j];
+			output[i] = raw; // raw vector to list
+		}
 	}
 	return output;
 }
@@ -725,6 +771,13 @@ Rcpp::List get_dim_sfc(Rcpp::List sfc) {
 			Rcpp::Named("_cls") = Rcpp::CharacterVector::create("XY"),
 			Rcpp::Named("_dim") = Rcpp::IntegerVector::create(2)
 		);
+	if (sfc.attr("points") != R_NilValue) {
+		Rcpp::NumericMatrix m = sfc.attr("points");
+		return Rcpp::List::create(
+			Rcpp::Named("_cls") = sfc.attr("points_dim"),
+			Rcpp::Named("_dim") = Rcpp::IntegerVector::create(m.ncol())
+		);
+	}
 
 	// we have data:
 	Rcpp::CharacterVector cls = sfc.attr("class");
diff --git a/tests/testthat/test-postgis_RPostgres.R b/tests/testthat/test-postgis_RPostgres.R
index b0e1ce4bf..feddc22b5 100644
--- a/tests/testthat/test-postgis_RPostgres.R
+++ b/tests/testthat/test-postgis_RPostgres.R
@@ -140,7 +140,7 @@ test_that("sf can write non-sf tables with geometries", {
     df <- as.data.frame(pts)
     expect_silent(st_write(df, pg, "df"))
     expect_silent(dfx <- st_read(pg, "df"))
-    expect_equal(df[["geometry"]], dfx[["geometry"]])
+    expect_equal(df$geometry, dfx$geometry)
     expect_silent(DBI::dbRemoveTable(pg, "df"))
 })
 
diff --git a/tests/testthat/test-sfc.R b/tests/testthat/test-sfc.R
index 5546d7930..0634792c0 100644
--- a/tests/testthat/test-sfc.R
+++ b/tests/testthat/test-sfc.R
@@ -120,10 +120,17 @@ test_that("st_is_longlat warns on invalid bounding box", {
 	expect_warning(st_is_longlat(st_sfc(st_point(c(0,-95)), crs = 4326)))
 })
 
+test_that("value replacement works for sfc_POINT",{
+	pts1<-st_geometry(st_as_sf(data.frame(x=1:3,y=1:3), coords = c("x","y")))
+	pts2<-st_geometry(st_as_sf(data.frame(x=4:5,y=4:5), coords = c("x","y")))
+	expect_identical(replace(pts1[],2:3,pts2), replace(pts1[],2:3,pts2[]))
+	expect_identical(replace(pts1,2:3,pts2[])[], replace(pts1[],2:3,pts2[]))
+	expect_identical(replace(pts1,2:3,pts2)[], replace(pts1[],2:3,pts2[]))
+	expect_identical(st_bbox(replace(pts1,2:3,pts2)),
+					 st_bbox(replace(pts1[],2:3,pts2[])))# check if bbox is correct without realization
+})
 test_that("bounding box is flipped when geometry is flipped", {
 	foo <- st_bbox(c(xmin = 0, xmax = 100, ymin = 0, ymax = 200)) |> st_as_sfc()
 	bar <- foo * matrix(c(1,0,0,-1), nrow = 2)
 	expect_equal(st_bbox(bar), st_bbox(c(xmin=0, ymin=-200, xmax=100, ymax=0)))
 })
-
-
diff --git a/tests/testthat/test-tidyverse-vctrs.R b/tests/testthat/test-tidyverse-vctrs.R
index 5c2fd8232..42c6dfa24 100644
--- a/tests/testthat/test-tidyverse-vctrs.R
+++ b/tests/testthat/test-tidyverse-vctrs.R
@@ -8,6 +8,12 @@ test_that("`sfc` vectors can be sliced", {
 	expect_identical(vctrs::vec_slice(x, 0), x[0])
 })
 
+test_that("`sfc` vectors with point matrix can be sliced", {
+	x = st_geometry(st_as_sf(data.frame(x=1:2, y=4:3), coords = 1:2))
+	expect_identical(vctrs::vec_slice(x, 1), x[1])
+	expect_equal(vctrs::vec_slice(x, 0), st_sfc(x[0], recompute_bbox=TRUE))
+})
+
 test_that("`n_empty` attribute of `sfc` vectors is restored", {
 	pt1 = st_sfc(st_point(c(NA_real_, NA_real_)))
 	pt2 = st_sfc(st_point(0:1))
diff --git a/tests/testthat/test-tm.R b/tests/testthat/test-tm.R
index 1a5e983c3..e73dc7a88 100644
--- a/tests/testthat/test-tm.R
+++ b/tests/testthat/test-tm.R
@@ -9,6 +9,7 @@ test_that("st_read and write handle date and time", {
 
   st_write(x[-4], shp[1], quiet = TRUE)
   x2 = st_read(shp[1], quiet = TRUE)
+  x2$geometry = x2$geometry[] # realize
   expect_equal(x[-4], x2, check.attributes=FALSE)
 # WKT2 CRS do not roundtrip for ESRI Shapefile
 
@@ -29,6 +30,7 @@ test_that("st_read and write handle date and time", {
 
   st_write(x[-4], shp[1], quiet = TRUE)
   x2 = st_read(shp[1], quiet = TRUE)
+  x2$geometry = x2$geometry[] # realize
   expect_equal(x[-4], x2, check.attributes=FALSE)
 # WKT2 CRS do not roundtrip for ESRI Shapefile
 
diff --git a/tests/testthat/test-wkb.R b/tests/testthat/test-wkb.R
index a5b10ca65..f11c2da20 100644
--- a/tests/testthat/test-wkb.R
+++ b/tests/testthat/test-wkb.R
@@ -34,7 +34,7 @@ test_that("Reading of big-endian and little-endian gives the same result", {
   y = structure(list("0x00200000010000714041061A800000000041145CAC00000000"), class = "WKB")
   expect_identical(st_as_sfc(x, EWKB = TRUE), st_as_sfc(y, EWKB = TRUE))
   expect_identical(st_as_sfc(x, EWKB = TRUE, pureR = TRUE), st_as_sfc(y, EWKB = TRUE, pureR = TRUE))
-  expect_identical(st_as_sfc(x, EWKB = TRUE), st_as_sfc(y, EWKB = TRUE, pureR = TRUE))
+  expect_identical(st_as_sfc(x, EWKB = TRUE)[], st_as_sfc(y, EWKB = TRUE, pureR = TRUE)) # [] realize
 })
 
 test_that("Reading of truncated buffers results in a proper error", {
diff --git a/tests/wkb.R b/tests/wkb.R
index 82f1e0b46..87e0bf76e 100644
--- a/tests/wkb.R
+++ b/tests/wkb.R
@@ -6,7 +6,7 @@ round_trip = function(x, EWKB = FALSE, pureR = FALSE) {
 	class(wkb) = "WKB"
 	# print(wkb)
 	y = st_as_sfc(wkb, EWKB = EWKB, pureR = pureR)
-	a = all.equal(x, y)
+	a = all.equal(x[], y[]) # realize both
 	if (length(a) == 1 && is.logical(a) && a)
 		TRUE
 	else {
diff --git a/tests/wkb.Rout.save b/tests/wkb.Rout.save
index 75e6dda4c..fd6738323 100644
--- a/tests/wkb.Rout.save
+++ b/tests/wkb.Rout.save
@@ -1,6 +1,6 @@
 
-R version 3.6.2 (2019-12-12) -- "Dark and Stormy Night"
-Copyright (C) 2019 The R Foundation for Statistical Computing
+R version 4.2.2 Patched (2022-11-10 r83330) -- "Innocent and Trusting"
+Copyright (C) 2022 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu (64-bit)
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -23,7 +23,7 @@ Type 'q()' to quit R.
 + 	class(wkb) = "WKB"
 + 	# print(wkb)
 + 	y = st_as_sfc(wkb, EWKB = EWKB, pureR = pureR)
-+ 	a = all.equal(x, y)
++ 	a = all.equal(x[], y[]) # realize both
 + 	if (length(a) == 1 && is.logical(a) && a)
 + 		TRUE
 + 	else {
@@ -99,4 +99,4 @@ POINT (0 0)
 > 
 > proc.time()
    user  system elapsed 
-  0.505   0.036   0.532 
+  1.022   0.741   0.874