R/plot.R

In edzer/sfr: Simple Features for R

kw_dflt = function(x, key.pos) {
	if (is.null(key.pos) || key.pos[1] == 0) # no key:
		return(lcm(0))

	font_scale = par("ps") / 12
	if (key.pos[1] == -1)
		lcm(1.8 * font_scale)
	else if (key.pos[1] %in% c(2, 4) && (is.character(x[[1]]) || is.factor(x[[1]]))) {
		strings = if (is.factor(x[[1]]))
				levels(x[[1]])
			else
				x[[1]]
		lcm(cm(max(strwidth(strings, "inches"))) * 1.3 + font_scale) # cm
		#lcm(cm(max(strwidth(strings, "inches"))) * 1.3) # cm
	} else
		lcm(1.8 * font_scale)
}


#' plot sf object
#'
#' plot one or more attributes of an sf object on a map
#' Plot sf object
#'
#' @param x object of class sf
#' @param y ignored
#' @param ... further specifications, see \link{plot_sf} and \link{plot} and details.
#' @param main title for plot (\code{NULL} to remove)
#' @param pal palette function, similar to \link{rainbow}, or palette values; if omitted, \code{sf.colors} is used
#' @param nbreaks number of colors breaks (ignored for \code{factor} or \code{character} variables)
#' @param breaks either a numeric vector with the actual breaks, or a name of a method accepted by the \code{style} argument of \link[classInt]{classIntervals}
#' @param max.plot integer; lower boundary to maximum number of attributes to plot; the default value (9) can be overridden by setting the global option \code{sf_max.plot}, e.g. \code{options(sf_max.plot=2)}
#' @param key.pos numeric; side to plot a color key: 1 bottom, 2 left, 3 top, 4 right; set to \code{NULL} to omit key completely, 0 to only not plot the key, or -1 to select automatically. If multiple columns are plotted in a single function call by default no key is plotted and every submap is stretched individually; if a key is requested (and \code{col} is missing) all maps are colored according to a single key. Auto select depends on plot size, map aspect, and, if set, parameter \code{asp}. If it has lenght 2, the second value, ranging from 0 to 1, determines where the key is placed in the available space (default: 0.5, center).
#' @param key.width amount of space reserved for the key (incl. labels), thickness/width of the scale bar
#' @param key.length amount of space reserved for the key along its axis, length of the scale bar
#' @param pch plotting symbol
#' @param cex symbol size
#' @param bg symbol background color
#' @param lty line type
#' @param lwd line width
#' @param col color for plotting features; if \code{length(col)} does not equal 1 or \code{nrow(x)}, a warning is emitted that colors will be recycled. Specifying \code{col} suppresses plotting the legend key.
#' @param border color of polygon border(s); using \code{NA} hides them
#' @param add logical; add to current plot? Note that when using \code{add=TRUE}, you may have to set \code{reset=FALSE} in the first plot command.
#' @param type plot type: 'p' for points, 'l' for lines, 'b' for both
#' @param reset logical; if \code{FALSE}, keep the plot in a mode that allows adding further map elements; if \code{TRUE} restore original mode after plotting \code{sf} objects with attributes; see details.
#' @param logz logical; if \code{TRUE}, use log10-scale for the attribute variable. In that case, \code{breaks} and \code{at} need to be given as log10-values; see examples.
#' @param extent object with an \code{st_bbox} method to define plot extent; defaults to \code{x}
#' @param xlim numeric; x-axis limits; overrides \code{extent}
#' @param ylim numeric; y-axis limits; overrides \code{extent}
#' @param compact logical; compact sub-plots over plotting space?
#' @method plot sf
#' @name plot
#' @details \code{plot.sf} maximally plots \code{max.plot} maps with colors following from attribute columns,
#' one map per attribute. It uses \code{sf.colors} for default colors. For more control over placement of individual maps,
#' set parameter \code{mfrow} with \link{par} prior to plotting, and plot single maps one by one; note that this only works
#' in combination with setting parameters \code{key.pos=NULL} (no legend) and \code{reset=FALSE}.
#'
#' \code{plot.sfc} plots the geometry, additional parameters can be passed on
#' to control color, lines or symbols.
#'
#' When setting \code{reset} to \code{FALSE}, the original device parameters are lost, and the device must be reset using \code{dev.off()} in order to reset it.
#'
#' parameter \code{at} can be set to specify where labels are placed along the key; see examples.
#' 
#' The features are plotted in the order as they apppear in the sf object. See examples for when a different plotting order is wanted.
#'
#' @examples
#' nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)
#' # plot single attribute, auto-legend:
#' plot(nc["SID74"])
#' # plot multiple:
#' plot(nc[c("SID74", "SID79")]) # better use ggplot2::geom_sf to facet and get a single legend!
#' # adding to a plot of an sf object only works when using reset=FALSE in the first plot:
#' plot(nc["SID74"], reset = FALSE)
#' plot(st_centroid(st_geometry(nc)), add = TRUE)
#' # log10 z-scale:
#' plot(nc["SID74"], logz = TRUE, breaks = c(0,.5,1,1.5,2), at = c(0,.5,1,1.5,2))
#' # and we need to reset the plotting device after that, e.g. by
#' layout(1)
#' # when plotting only geometries, the reset=FALSE is not needed:
#' plot(st_geometry(nc))
#' plot(st_geometry(nc)[1], col = 'red', add = TRUE)
#' # add a custom legend to an arbitray plot:
#' layout(matrix(1:2, ncol = 2), widths = c(1, lcm(2)))
#' plot(1)
#' .image_scale(1:10, col = sf.colors(9), key.length = lcm(8), key.pos = 4, at = 1:10)
#' # manipulate plotting order, plot largest polygons first:
#' p = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
#' x = st_sf(a=1:4, st_sfc(p, p * 2, p * 3, p * 4)) # plot(x, col=2:5) only shows the largest polygon!
#' plot(x[order(st_area(x), decreasing = TRUE),], col = 2:5) # plot largest polygons first
#' 
#' @export
plot.sf <- function(x, y, ..., main, pal = NULL, nbreaks = 10, breaks = "pretty",
		max.plot = getOption("sf_max.plot", default = 9),
		key.pos = get_key_pos(x, ...), key.length = .618, key.width = kw_dflt(x, key.pos),
		reset = TRUE, logz = FALSE, extent = x, xlim = st_bbox(extent)[c(1,3)],
		ylim = st_bbox(extent)[c(2,4)], compact = FALSE) {

	stopifnot(missing(y))
	nbreaks.missing = missing(nbreaks)
	key.pos.missing = missing(key.pos)
	max_plot_missing = missing(max.plot)
	dots = list(...)
	col_missing = is.null(dots$col)
	breaks_numeric = is.numeric(breaks)
	reset_layout_needed = reset

	x = swap_axes_if_needed(x)

	opar = par(no.readonly = TRUE)
	if (ncol(x) > 2 && !isTRUE(dots$add)) { # multiple maps to plot...
		cols = setdiff(names(x), attr(x, "sf_column"))
		lt = .get_layout(st_bbox(x), min(max.plot, length(cols)), par("din"), key.pos[1], key.width)
		if (key.pos.missing || key.pos[1] == -1)
			key.pos = lt$key.pos
		layout(lt$m, widths = lt$widths, heights = lt$heights, respect = compact)

		if (isTRUE(dots$axes))
			par(mar = c(2.1, 2.1, 1.2, 0))
		else
			par(mar = c(0, 0, 1.2, 0))

		if (max_plot_missing)
			max.plot = prod(lt$mfrow)

		if (isTRUE(is.finite(max.plot)) && ncol(x) - 1 > max.plot &&
				max_plot_missing && is.null(options("sf_max.plot")[[1]]))
			warning(paste("plotting the first", max.plot, "out of", ncol(x)-1,
				"attributes; use max.plot =", ncol(x) - 1, "to plot all"), call. = FALSE)

		# col selection may have changed; set cols again:
		cols = setdiff(names(x), attr(x, "sf_column"))
		if (length(cols) > max.plot)
			cols = cols[1:max.plot]

		if (!is.null(key.pos)) {
			values = do.call(c, as.data.frame(x)[cols])
			if (is.character(values))
				values = as.factor(values)
			if (logz)
				values = log10(values)
			if (is.character(breaks) && is.numeric(values)) { # compute breaks from values:
				v0 = values[!is.na(values)]
				n.unq = length(unique(v0))
				breaks = if (! all(is.na(values)) && n.unq > 1)
						classInt::classIntervals(v0, min(nbreaks, n.unq),
							breaks, warnSmallN = FALSE)$brks
					else
						range(values, na.rm = TRUE) # lowest and highest!
				nbreaks = length(breaks) - 1
			}
		}
		if (nbreaks.missing && is.numeric(breaks))
			nbreaks = length(breaks) - 1

		# loop over each map to plot:
		lapply(cols, function(cname) plot(x[, cname], main = cname,
			pal = pal, nbreaks = nbreaks, breaks = breaks, key.pos = NULL, reset = FALSE,
			logz = logz, xlim = xlim, ylim = ylim,...))

		for (i in seq_len(prod(lt$mfrow) - length(cols))) # empty panels:
			plot.new()

		# plot key?
		if (!is.null(key.pos) && key.pos[1] != 0 && col_missing) {
			if (is.null(pal))
				pal = function(n) sf.colors(n, categorical = is.factor(values))
			colors = if (is.function(pal))
					pal(nbreaks)
				else
					pal
			if (is.factor(values))
				.image_scale_factor(levels(values), colors, key.pos = key.pos,
					key.width = key.width, key.length = key.length, ...)
			else
				.image_scale(values, colors, breaks = breaks, key.pos = key.pos,
					key.length = key.length, logz = logz, ...)
		}

	} else { # single map, or dots$add == TRUE:
		if (ncol(x) == 1) { # no attributes to choose colors from: plot geometry
			plot(st_geometry(x), xlim = xlim, ylim = ylim, ...)
			reset_layout_needed = FALSE
		} else { # generate plot with colors and possibly key
			if (ncol(x) > 2) { # add = TRUE
				warning("ignoring all but the first attribute")
				x = x[,1]
			}
			# store attribute in "values":
			values = x[[setdiff(names(x), attr(x, "sf_column"))]]

			if (is.list(values))
				stop("plotting list-columns not supported") # nocov

			if (is.character(values))
				values = as.factor(values)
			else if (logz)
				values = log10(as.numeric(values))

			if (is.null(pal))
				pal = function(n) sf.colors(n, categorical = is.factor(values))
			else if (! col_missing)
				stop("specify only one of `col' and `pal'")

			if (col_missing) { # compute colors from values:
				col = if (is.factor(values)) {
						if (key.pos.missing && nlevels(values) > 30) # doesn't make sense:
							key.pos = NULL
						colors = if (is.function(pal))
								pal(nlevels(values))
							else
								pal
						colors[as.numeric(values)]
					} else {
						if (! inherits(values, c("POSIXt", "Date")))
							values = as.numeric(values) # drop units, if any
						if (is.character(breaks)) { # compute breaks from values:
							v0 = values[!is.na(values)]
							n.unq = length(unique(v0))
							breaks = if (! all(is.na(values)) && n.unq > 1)
									classInt::classIntervals(v0, min(nbreaks, n.unq),
										breaks, warnSmallN = FALSE)$brks
								else
									range(values, na.rm = TRUE) # lowest and highest!
						}
						# this is necessary if breaks were specified either as character or as numeric
						# "pretty" takes nbreaks as advice only:
						nbreaks = length(breaks) - 1

						cuts = if (all(is.na(values)))
								rep(NA_integer_, length(values))
							else if (!breaks_numeric && diff(range(values, na.rm = TRUE)) == 0)
								ifelse(is.na(values), NA_integer_, 1L)
							else
								cut(values, breaks, include.lowest = TRUE)
						colors = if (is.function(pal))
								pal(nbreaks)
							else
								pal
						colors[cuts]
					}
			} else {
				col = dots$col
				if (length(col) != 1 && length(col) != nrow(x))
					warning("col is not of length 1 or nrow(x): colors will be recycled; use pal to specify a color palette")
				key.pos = NULL # no key!
			}

			if (!isTRUE(dots$add) && !is.null(key.pos) && !all(is.na(values)) &&
					(is.factor(values) || length(unique(na.omit(values))) > 1 || breaks_numeric) && # 2065
					length(col) > 1) { # plot key?

				switch(key.pos[1],
					layout(matrix(c(2,1), nrow = 2, ncol = 1),
						widths = 1, heights = c(1, key.width)), # 1 bottom
					layout(matrix(c(1,2), nrow = 1, ncol = 2),
						widths = c(key.width, 1), heights = 1), # 2 left
					layout(matrix(c(1,2), nrow = 2, ncol = 1),
						widths = 1, heights = c(key.width, 1)), # 3 top
					layout(matrix(c(2,1), nrow = 1, ncol = 2),
						widths = c(1, key.width), heights = 1)  # 4 right
				)

				if (is.factor(values)) {
					.image_scale_factor(levels(values), colors, key.pos = key.pos,
						key.width = key.width, key.length = key.length, ...)
				} else
					.image_scale(values, colors, breaks = breaks, key.pos = key.pos,
						key.length = key.length, logz = logz, ...)
			} else
				reset_layout_needed = FALSE # as we didn't call layout()
			# plot the map:
			if (!isTRUE(dots$add)) {
				mar = c(1, 1, 1.2, 1)
				if (isTRUE(dots$axes))
					mar[1:2] = 2.1
				par(mar = mar)
			}
			if (col_missing)
				plot(st_geometry(x), col = col, xlim = xlim, ylim = ylim, ...)
			else
				plot(st_geometry(x), xlim = xlim, ylim = ylim, ...)
		}
		if (! isTRUE(dots$add)) { # title?
			if (missing(main)) {
				main = setdiff(names(x), attr(x, "sf_column"))
				if (length(main) && inherits(x[[main]], "units"))
					main = make_unit_label(main, x[[main]])
			}
			localTitle <- function(..., extent, col, bg, pch, cex, lty, lwd, axes, type, bgMap,
					border, graticule, xlim, ylim, asp, bgc, xaxs, yaxs, lab, setParUsrBB,
					expandBB, col_graticule, at, lon, lat, crs, datum, ndiscr, margin) # absorb
				title(...)
			localTitle(main, ...)
		}
	}
	if (!isTRUE(dots$add) && reset) { # reset device: 
		if (reset_layout_needed) 
			layout(matrix(1))
		par(opar)
	} 
	invisible()
}

swap_axes_if_needed = function(x) {
	crs = st_crs(x)
	if (st_axis_order() && !is.na(crs) && crs$yx)
		st_transform(x, pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1")
	else
		x
}


#' @name plot
#' @export
get_key_pos = function(x, ...) {
	bb = st_bbox(x)
	if (anyNA(bb) || (inherits(x, "sf") && ncol(x) > 2))
		NULL
	else {
		pin = par("pin") # (width, height)
		asp_plt = pin[2]/pin[1] # y/x: < 1 means wide
		asp_box = diff(bb[c(4,2)]) / diff(bb[c(3,1)])
		asp = list(...)$asp
		if (is.null(asp))
			asp <- ifelse(isTRUE(st_is_longlat(x)), 1/cos((mean(bb[c(2,4)]) * pi)/180), 1.0)
		asp_box = asp_box * asp
		if (!is.finite(asp_box) || asp_box < asp_plt) # plot is wider than device: below
			1
		else # plot is taller than device: to the right
			4
	}
}

#' @name plot
#' @method plot sfc_POINT
#' @export
plot.sfc_POINT = function(x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1,
		type = 'p', add = FALSE) {
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	npts = length(x)
	pch = rep(pch, length.out = npts)
	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
	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],
			cex = cex[ne], lwd = lwd, lty = lty, type = type)
	}
}

#' @name plot
#' @method plot sfc_MULTIPOINT
#' @export
plot.sfc_MULTIPOINT = function(x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1,
		type = 'p', add = FALSE) {
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	n = length(x)
	pch = rep(pch, length.out = n)
	col = rep(col, length.out = n)
	bg = rep(bg, length.out = n)
	cex = rep(cex, length.out = n)
	lwd = rep(lwd, length.out = n)
	lty = rep(lty, length.out = n)
	non_empty = ! st_is_empty(x)
	lapply(seq_along(x), function(i)
	  if (non_empty[i])
		points(x[[i]], pch = pch[i], col = col[i], bg = bg[i],
			cex = cex[i], lwd = lwd[i], lty = lty[i], type = type))
	invisible(NULL)
}

#' @name plot
#' @method plot sfc_LINESTRING
#' @export
plot.sfc_LINESTRING = function(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = 'l',
		add = FALSE) {
# FIXME: take care of lend, ljoin, and lmitre
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	lty = rep(lty, length.out = length(x))
	lwd = rep(lwd, length.out = length(x))
	col = rep(col, length.out = length(x))
	pch  = rep(pch, length.out = length(x))
	non_empty = ! st_is_empty(x)
	lapply(seq_along(x), function(i)
	  if (non_empty[i])
		lines(x[[i]], lty = lty[i], lwd = lwd[i], col = col[i], pch = pch[i], type = type))
	invisible(NULL)
}

#' @name plot
#' @method plot sfc_CIRCULARSTRING
#' @export
plot.sfc_CIRCULARSTRING = function(x, y, ...) {
	plot(st_cast(x, "LINESTRING"),  ...)
}

#' @name plot
#' @method plot sfc_MULTILINESTRING
#' @export
plot.sfc_MULTILINESTRING = function(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = 'l',
		add = FALSE) {
# FIXME: take care of lend, ljoin, and lmitre
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	lty = rep(lty, length.out = length(x))
	lwd = rep(lwd, length.out = length(x))
	col = rep(col, length.out = length(x))
	pch  = rep(pch, length.out = length(x))
	non_empty = ! st_is_empty(x)
	lapply(seq_along(x), function(i)
	  if (non_empty[i])
		lapply(x[[i]], function(L)
			lines(L, lty = lty[i], lwd = lwd[i], col = col[i], pch = pch[i], type = type)))
	invisible(NULL)
}

# sf (list) -> polypath (mtrx) : rbind polygon rings with NA rows inbetween
p_bind = function(lst) {
	if (length(lst) == 1)
		lst[[1]]
	else {
		ret = vector("list", length(lst) * 2 - 1)
		ret[seq(1, length(lst) * 2 - 1, by = 2)] = lst # odd elements
		ret[seq(2, length(lst) * 2 - 1, by = 2)] = NA  # even elements
		do.call(rbind, ret) # replicates the NA to form an NA row
	}
}

#' @name plot
#' @param rule see \link[graphics]{polypath}; for \code{winding}, exterior ring direction should be opposite that of the holes; with \code{evenodd}, plotting is robust against misspecified ring directions
#' @param xpd see \link[graphics]{par}; sets polygon clipping strategy; only implemented for POLYGON and MULTIPOLYGON
#' @export
plot.sfc_POLYGON = function(x, y, ..., lty = 1, lwd = 1, col = NA, cex = 1, pch = NA, border = 1,
		add = FALSE, rule = "evenodd", xpd = par("xpd")) {
# FIXME: take care of lend, ljoin, and lmitre
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	lty = rep(lty, length.out = length(x))
	lwd = rep(lwd, length.out = length(x))
	col = rep(col, length.out = length(x))
	border = rep(border, length.out = length(x))
	non_empty = ! st_is_empty(x)
	lapply(seq_along(x), function(i)
	  if (non_empty[i])
		polypath(p_bind(x[[i]]), border = border[i], lty = lty[i], lwd = lwd[i], col = col[i], rule = rule, xpd = xpd))
#	if (any(!is.na(pch))) {
#		pch = rep(pch, length.out = length(x))
#		cex = rep(cex, length.out = length(x))
#		lapply(seq_along(x), function(i)
#		  if (non_empty[i])
#			points(p_bind(x[[i]]), pch = pch[i], cex = cex[i], type = 'p'))
#	}
	invisible(NULL)
}

#' @name plot
#' @method plot sfc_MULTIPOLYGON
#' @export
plot.sfc_MULTIPOLYGON = function(x, y, ..., lty = 1, lwd = 1, col = NA, border = 1, add = FALSE, 
		rule = "evenodd", xpd = par("xpd")) {
# FIXME: take care of lend, ljoin, and lmitre
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	lty = rep(lty, length.out = length(x))
	lwd = rep(lwd, length.out = length(x))
	col = rep(col, length.out = length(x))
	border = rep(border, length.out = length(x))
	non_empty = ! st_is_empty(x)
	lapply(seq_along(x), function(i)
	  if (non_empty[i])
		lapply(x[[i]], function(L)
			polypath(p_bind(L), border = border[i], lty = lty[i], lwd = lwd[i], col = col[i], rule = rule, xpd = xpd)))
	invisible(NULL)
}

# plot single geometrycollection:
plot_gc = function(x, pch, cex, bg, border = 1, lty, lwd, col, add) {
	lapply(x, function(subx) {
		args = list(st_sfc(subx), pch = pch, cex = cex, bg = bg, border = border,
			lty = lty, lwd = lwd, col = col, add = TRUE)
		fn = switch(class(subx)[2],
			POINT = plot.sfc_POINT,
			MULTIPOINT = plot.sfc_MULTIPOINT,
			LINESTRING = plot.sfc_LINESTRING,
			MULTILINESTRING = plot.sfc_MULTILINESTRING,
			POLYGON = plot.sfc_POLYGON,
			CIRCULARSTRING = plot.sfc_CIRCULARSTRING,
			MULTIPOLYGON = plot.sfc_MULTIPOLYGON,
			MULTISURFACE = plot.sfc_GEOMETRYCOLLECTION,
			CURVEPOLYGON = plot.sfc_GEOMETRYCOLLECTION,
			COMPOUNDCURVE = plot.sfc_GEOMETRYCOLLECTION,
			GEOMETRYCOLLECTION = plot.sfc_GEOMETRYCOLLECTION,
			stop(paste("plotting of", class(x)[2], "not yet supported: use st_cast?"))
		)
		do.call(fn, args)
	})
	invisible(NULL)
}

#' @name plot
#' @method plot sfc_GEOMETRYCOLLECTION
#' @export
plot.sfc_GEOMETRYCOLLECTION = function(x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1,
	col = 1, border = 1, add = FALSE) {
# FIXME: take care of lend, ljoin, xpd, and lmitre
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	cex = rep(cex, length.out = length(x))
	pch = rep(pch, length.out = length(x))
	lty = rep(lty, length.out = length(x))
	lwd = rep(lwd, length.out = length(x))
	col = rep(col, length.out = length(x))
	border = rep(border, length.out = length(x))
	lapply(seq_along(x), function(i) plot_gc(x[[i]],
			pch = pch[i], cex = cex[i], bg = bg[i], border = border[i], lty = lty[i],
			lwd = lwd[i], col = col[i]))
	invisible(NULL)
}

#' @name plot
#' @method plot sfc_GEOMETRY
#' @export
plot.sfc_GEOMETRY = function(x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1,
	col = ifelse(st_dimension(x) == 2, NA, 1), border = 1, add = FALSE) {
	stopifnot(missing(y))
	if (! add)
		plot_sf(x, ...)
	cex = rep(cex, length.out = length(x))
	pch = rep(pch, length.out = length(x))
	lty = rep(lty, length.out = length(x))
	lwd = rep(lwd, length.out = length(x))
	col = rep(col, length.out = length(x))
	border = rep(border, length.out = length(x))
	lapply(seq_along(x), function(i) plot_gc(st_sfc(x[[i]]),
			pch = pch[i], cex = cex[i], bg = bg[i], border = border[i], lty = lty[i],
			lwd = lwd[i], col = col[i]))
	invisible(NULL)
}

#' @name plot
#' @method plot sfg
#' @export
plot.sfg = function(x, ...) {
	plot(st_sfc(x), ...)
}

# set up plotting area & axes; reuses sp:::plot.Spatial
#' @name plot
#' @param xlim see \link{plot.window}
#' @param ylim see \link{plot.window}
#' @param asp see below, and see \link{par}
#' @param axes logical; should axes be plotted? (default FALSE)
#' @param bgc background color
#' @param xaxs see \link{par}
#' @param yaxs see \link{par}
#' @param lab see \link{par}
#' @param setParUsrBB default FALSE; set the \code{par} \dQuote{usr} bounding box; see below
#' @param bgMap object of class \code{ggmap}, or returned by function \code{RgoogleMaps::GetMap}
#' @param expandBB numeric; fractional values to expand the bounding box with,
#' in each direction (bottom, left, top, right)
#' @param graticule logical, or object of class \code{crs} (e.g., \code{st_crs(4326)} for a WGS84 graticule), or object created by \link{st_graticule}; \code{TRUE} will give the WGS84 graticule
#' or object returned by \link{st_graticule}
#' @param col_graticule color to used for the graticule (if present)
#' @export
#' @details \code{plot_sf} sets up the plotting area, axes, graticule, or webmap background; it
#' is called by all \code{plot} methods before anything is drawn.
#'
#' The argument \code{setParUsrBB} may be used to pass the logical value \code{TRUE} to functions within \code{plot.Spatial}. When set to \code{TRUE}, par(\dQuote{usr}) will be overwritten with \code{c(xlim, ylim)}, which defaults to the bounding box of the spatial object. This is only needed in the particular context of graphic output to a specified device with given width and height, to be matched to the spatial object, when using par(\dQuote{xaxs}) and par(\dQuote{yaxs}) in addition to \code{par(mar=c(0,0,0,0))}.
#'
#' The default aspect for map plots is 1; if however data are not
#' projected (coordinates are long/lat), the aspect is by default set to
#' 1/cos(My * pi/180) with My the y coordinate of the middle of the map
#' (the mean of \code{ylim}, which defaults to the y range of bounding box). This
#' implies an \href{https://en.wikipedia.org/wiki/Equirectangular_projection}{Equirectangular projection}.
#'
plot_sf = function(x, xlim = NULL, ylim = NULL, asp = NA, axes = FALSE, bgc = par("bg"), ...,
    xaxs, yaxs, lab, setParUsrBB = FALSE, bgMap = NULL, expandBB = c(0,0,0,0), graticule = NA_crs_,
	col_graticule = 'grey', border, extent = x) {

# sp's bbox: matrix
#   min max
# x
# y

	bbox = matrix(st_bbox(extent), 2, dimnames = list(c("x", "y"), c("min", "max")))
	# expandBB: 1=below, 2=left, 3=above and 4=right.
	expBB = function(lim, expand) c(lim[1] - expand[1] * diff(lim), lim[2] + expand[2] * diff(lim))
	if (is.null(xlim))
		xlim <- expBB(bbox[1,], expandBB[c(2,4)])
	if (is.null(ylim))
		ylim <- expBB(bbox[2,], expandBB[c(1,3)])
	if (is.na(asp))
		asp <- ifelse(isTRUE(st_is_longlat(x)), 1/cos((mean(ylim) * pi)/180), 1.0)

	if (anyNA(bbox))
		stop("NA value(s) in bounding box. Trying to plot empty geometries?")

	plot.new()

	args = list(xlim = xlim, ylim = ylim, asp = asp)
	if (!missing(xaxs)) args$xaxs = xaxs
	if (!missing(yaxs)) args$yaxs = yaxs
	if (!missing(lab)) args$lab = lab
	do.call(plot.window, args)

	if (setParUsrBB)
		par(usr=c(xlim, ylim))
	pl_reg <- par("usr")
	rect(xleft = pl_reg[1], ybottom = pl_reg[3], xright = pl_reg[2],
		ytop = pl_reg[4], col = bgc, border = FALSE)
	linAxis = function(side, ..., lon, lat, ndiscr, reset, at) axis(side = side, ...)
	if (! missing(graticule)) {
		g = if (isTRUE(graticule))
				st_graticule(pl_reg[c(1,3,2,4)], st_crs(x), st_crs(4326), ...)
			else if (inherits(graticule, "crs") && !is.na(graticule))
				st_graticule(pl_reg[c(1,3,2,4)], st_crs(x), graticule, ...)
			else
				graticule
		plot(st_geometry(g), col = col_graticule, add = TRUE)
		box()
		if (axes) {
			sel = g$type == "E" & g$plot12
			linAxis(1L, g$x_start[sel], parse(text = g$degree_label[sel]), ...)
			sel = g$type == "N" & g$plot12
			linAxis(2L, g$y_start[sel], parse(text = g$degree_label[sel]), ...)
		}
	} else if (axes) {
		box()
		if (isTRUE(st_is_longlat(x))) {
			local_degAxis = function(side, ..., at) .degAxis(side, ...) # absorb at
			local_degAxis(1, ...)
			local_degAxis(2, ...)
		} else {
			linAxis(1, ...)
			linAxis(2, ...)
		}
	}
	localTitle <- function(..., col, bgc, pch, cex, lty, lwd, lon, lat, ndiscr, at, labels, reset) title(...)
	localTitle(...)
	if (!is.null(bgMap)) {
		mercator = FALSE
		if (inherits(bgMap, "ggmap")) {
			bb = bb2merc(bgMap, "ggmap")
			mercator = TRUE
		} else if (all(c("lat.center","lon.center","zoom","myTile","BBOX") %in% names(bgMap))) {
			# an object returned by RgoogleMaps::GetMap
			bb = bb2merc(bgMap, "RgoogleMaps")
			bgMap = bgMap$myTile
			mercator = TRUE
		} else
			bb = c(xlim[1], ylim[1], xlim[2], ylim[2]) # can be any crs!
		if (mercator &&  st_crs(x) != st_crs(3857))
			warning("crs of plotting object differs from that of bgMap, which is assumed to be st_crs(3857)") # nocov
		rasterImage(bgMap, bb[1], bb[2], bb[3], bb[4], interpolate = FALSE)
	}
}


#' @param n integer; number of colors
#' @param cutoff.tails numeric, in `[0,0.5]` start and end values
#' @param alpha numeric, in `[0,1]`, transparency
#' @param categorical logical; do we want colors for a categorical variable? (see details)
#' @name plot
#' @export
#' @details non-categorical colors from \code{sf.colors} were taken from \link[sp]{bpy.colors}, with modified \code{cutoff.tails} defaults
#' If categorical is \code{TRUE}, default colors are from \url{https://colorbrewer2.org/} (if n < 9, Set2, else Set3).
#' @examples
#' sf.colors(10)
sf.colors = function (n = 10, cutoff.tails = c(0.35, 0.2), alpha = 1, categorical = FALSE) {
	if (categorical) {
		cb = if (n <= 8)
		# 8-class Set2:
		c('#66c2a5','#fc8d62','#8da0cb','#e78ac3','#a6d854','#ffd92f','#e5c494','#b3b3b3')
		# 12-class Set3:
		else c('#8dd3c7','#ffffb3','#bebada','#fb8072','#80b1d3','#fdb462','#b3de69','#fccde5','#d9d9d9','#bc80bd','#ccebc5','#ffed6f')
		# TODO: deal with alpha
		if (alpha != 1.0)
			cb = paste0(cb, as.hexmode(ceiling(alpha * 255)))
		rep(cb, length.out = n)
	} else {
		i = seq(0.5 * cutoff.tails[1], 1 - 0.5 * cutoff.tails[2], length.out = n)
   		r = ifelse(i < .25, 0, ifelse(i < .57, i / .32 - .78125, 1))
   		g = ifelse(i < .42, 0, ifelse(i < .92, 2 * i - .84, 1))
   		b = ifelse(i < .25, 4 * i, ifelse(i < .42, 1,
			ifelse(i < .92, -2 * i + 1.84, i / .08 - 11.5)))
		rgb(r, g, b, alpha)
	}
}

# get the aspect ratio of a bounding box, for geodetic coords true scale at mid latitude:
get_asp = function(bb) {
	asp = diff(bb[c(2,4)])/diff(bb[c(1,3)])
	if (!is.finite(asp)) # 0/0
		asp = 1
	if (isTRUE(st_is_longlat(bb)))
		asp = asp / cos(mean(bb[c(2,4)]) * pi /180)
	asp
}

#' functions only exported to be used internally by stars
#' @keywords internal
#' @name stars
#' @param bb ignore
#' @param n ignore
#' @param total_size ignore
#' @param key.width ignore
#' @param key.length ignore
#' @param mfrow length-2 integer vector with number of rows, columns
#' @param main main or sub title
#' @export
.get_layout = function(bb, n, total_size, key.pos, key.width, mfrow = NULL, main = NULL) {
# return list with "m" matrix, "key.pos", "widths" and "heights" fields
# if key.pos = -1 on input, it will be a return value, "optimally" placed
	asp = get_asp(bb)
	strip = if (is.character(main))
			# strheight(main, "inches") 
			par("cin")[2]
		else 
			0.0
	size = function(nrow, n, asp, strip = 0) { # given nrow n asp, what size does a single tile occupy?
		ncol = ceiling(n / nrow)
		xsize = total_size[1] / ncol
		ysize = xsize  * asp + strip
		if (xsize * ysize * n > prod(total_size)) {
			ysize = total_size[2] / nrow - strip
			xsize = ysize / asp
		}
		c(xsize, ysize)
	}
	sz = vapply(1:n, function(nrow) size(nrow, n, asp, strip), c(0.0, 0.0))
	if (is.null(mfrow)) {
		nrow = which.max(apply(sz, 2, prod))
		ncol = ceiling(n / nrow)
	} else {
		stopifnot(is.numeric(mfrow), length(mfrow) == 2)
		nrow = mfrow[1]
		ncol = mfrow[2]
	}
	xsize = sz[1, nrow]
	ysize = sz[2, nrow]
	asp = ysize / xsize

	ret = list()
	ret$mfrow = c(nrow, ncol)

	# the following is right now only used by stars; FIXME:
	# nocov start
	ret$key.pos = if (!is.null(key.pos) && key.pos == -1L) { # figure out here: right or bottom?
			newasp = asp * nrow / ncol # of the composition
			dispasp = total_size[2] / total_size[1]
			ifelse(newasp > dispasp, 4, 1)
		} else
			key.pos

	m = matrix(seq_len(nrow * ncol), nrow, ncol, byrow = TRUE)
	if (!is.null(ret$key.pos) && ret$key.pos != 0) { # add key row or column:
		k = key.width
		n = nrow * ncol + 1
		switch(ret$key.pos,
			{ ret$m = rbind(m, n); ret$widths = c(rep(1, ncol)); ret$heights = c(rep(asp, nrow), k) },
			{ ret$m = cbind(n, m); ret$widths = c(k, rep(1, ncol)); ret$heights = c(rep(asp, nrow)) },
			{ ret$m = rbind(n, m); ret$widths = c(rep(1, ncol)); ret$heights = c(k, rep(asp, nrow)) },
			{ ret$m = cbind(m, n); ret$widths = c(rep(1, ncol), k); ret$heights = c(rep(asp, nrow)) }
		)
	} else {
		ret$m = m
		ret$widths = rep(1, ncol)
		ret$heights = rep(asp, nrow)
	}
	# nocov end
	ret
}

bb2merc = function(x, cls = "ggmap") { # return bbox in the appropriate "web mercator" CRS
	wgs84 = st_crs(4326)
	merc =  st_crs(3857) # http://wiki.openstreetmap.org/wiki/EPSG:3857
	pts = if (cls == "ggmap") {
		b = vapply(attr(x, "bb"), c, 0.0)
		st_sfc(st_point(c(b[2:1])), st_point(c(b[4:3])), crs = wgs84)
	} else if (cls == "RgoogleMaps")
		st_sfc(st_point(rev(x$BBOX$ll)), st_point(rev(x$BBOX$ur)), crs = wgs84)
	else
		stop("unknown cls")
	st_bbox(st_transform(pts, merc))
}

#' @rdname stars
#' @param side ignore
#' @param at ignore
#' @param labels ignore
#' @param lon ignore
#' @param lat ignore
#' @param ndiscr ignore
#' @param reset ignore
#' @export
.degAxis = function (side, at, labels, ..., lon, lat, ndiscr, reset) {
	if (missing(at))
       	at = axTicks(side)
	if (missing(labels)) {
		labels = FALSE
		if (side == 1 || side == 3)
			labels = parse(text = degreeLabelsEW(at))
		else if (side == 2 || side == 4)
			labels = parse(text = degreeLabelsNS(at))
	}
	axis(side, at = at, labels = labels, ...)
}

# find out where to place the legend key:
# given range r = (a, b), key.length l, key offset o, return a value range that:
#  (i)  scales such that (b - a) / (y - x) = l, and
#  (ii) shifts linearly within [x, y] from a = x when o = 0 to b = y when o = 1
xy_from_r = function(r, l, o) {
	stopifnot(length(r) == 2, l <= 1, l > 0, o >= 0, o <= 1)
	r = as.numeric(r)
	a = r[1]; b = r[2]
	if (o == 1) {
		y = b
		x = b - (b - a)/l
	} else {
		i = o / (o - 1)
		y = (a + (b - a)/l - i * b)/(1 - i)
		x = i * (y - b) + a
	}
	c(x, y)
}

#' @rdname stars
#' @param z ignore
#' @param col ignore
#' @param breaks ignore
#' @param key.pos ignore
#' @param add.axis ignore
#' @param axes ignore
#' @param logz ignore
#' @param ... ignore
#' @param lab ignore
#' @param cex.axis see \link{par}
#' @export
.image_scale = function(z, col, breaks = NULL, key.pos, add.axis = TRUE,
		at = NULL, ..., axes = FALSE, key.length, logz = FALSE, lab = "", 
		cex.axis = par("cex.axis")) {
	if (!is.null(breaks) && length(breaks) != (length(col) + 1))
		stop("must have one more break than colour")
	stopifnot(is.null(lab) || is.character(lab) || is.expression(lab))
	lab_set = (is.character(lab) && lab != "") || is.expression(lab)
	zlim = range(z, na.rm = TRUE)
	if (is.null(breaks))
		breaks = seq(zlim[1], zlim[2], length.out = length(col) + 1)
	offset = 0.5
	if (length(key.pos) == 2) {
		offset = key.pos[2]
		key.pos = key.pos[1]
	}
	if (is.character(key.length)) {
		kl = as.numeric(gsub(" cm", "", key.length))
		sz = if (key.pos %in% c(1,3))
				dev.size("cm")[1]
			else
				dev.size("cm")[2]
		key.length = kl/sz
	}
	if (is.null(at)) {
		br = range(breaks)
		at = pretty(br)
		at = at[at > br[1] & at < br[2]]
	}
	if (key.pos %in% c(1,3)) {
		ylim = c(0, 1)
		xlim = xy_from_r(range(breaks), key.length, offset)
		mar = c(0, ifelse(axes, 2.1, 1), 0, 1)
	}
	if (key.pos %in% c(2,4)) {
		ylim = xy_from_r(range(breaks), key.length, offset)
		xlim = c(0, 1)
		mar = c(ifelse(axes, 2.1, 1), 0, 1.2, 0)
	}
	mar[key.pos] = 2.1 + 1.5 * lab_set
	par(mar = mar)

	plot(1, 1, type = "n", ylim = ylim, xlim = xlim, axes = FALSE,
		xlab = "", ylab = "", xaxs = "i", yaxs = "i")
	if (!is.null(lab) && lab != "")
		mtext(lab, side = key.pos, line = 2.5, cex = .8)
	poly = vector(mode="list", length(col))
	for (i in seq(poly))
		poly[[i]] = c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
	offset = 0.2
	offs = switch(key.pos,
		c(0, 0, -offset, -offset),
		c(0, 0, -offset, -offset),
		c(offset, offset, 0, 0),
		c(offset, offset, 0, 0))
	for(i in seq_along(poly)) {
		if (key.pos %in% c(1,3))
			polygon(poly[[i]], c(0, 0, 1, 1) + offs, col = col[i], border = NA)
		if (key.pos %in% c(2,4))
			polygon(c(0, 0, 1, 1) + offs, poly[[i]], col = col[i], border = NA)
	}

	# box() now would draw around [0,1]:
	bx = c(breaks[1], rep(tail(breaks, 1), 2), breaks[1])
	if (key.pos %in% c(1, 3))
		polygon(bx, c(0, 0, 1, 1) + offs, col = NA, border = 'black')
	if (key.pos %in% c(2, 4))
		polygon(c(0, 0, 1, 1) + offs, bx, col = NA, border = 'black')

	labels = if (logz)
			parse(text = paste0("10^", at))
		else if (inherits(breaks, c("POSIXt", "Date")))
			format(at)
		else
			TRUE

	if (add.axis)
		axis(key.pos, at = at, labels = labels, cex.axis = cex.axis)
}

#' @rdname stars
#' @param key.width ignore
#' @export
.image_scale_factor = function(z, col, key.pos, add.axis = TRUE,
	..., axes = FALSE, key.width, key.length, cex.axis = par("cex.axis")) {

	to_num = function(x) as.numeric(gsub(" cm", "", x))
	n = length(z)
	if ((kw <- to_num(key.width)) < to_num(kw_dflt(list(z), key.pos))) { # cut z until it fits:
		m = max(nchar(z))
		while(kw < to_num(kw_dflt(list(z), key.pos))) { # cut z:
			m = m - 1
			z = substr(z, 1, m)
		}
	}
	ksz = max(1.5 + max(nchar(z))/2, max(strwidth(z, "inches")) / par("cin")[1]) # in "mar" lines
	breaks = (0:n) + 0.5
	offset = 0.5
	if (length(key.pos) == 2) {
		offset = key.pos[2]
		key.pos = key.pos[1]
	}
	if (is.character(key.length)) {
		kl = to_num(key.length)
		sz = if (key.pos %in% c(1, 3))
				dev.size("cm")[1]
			else
				dev.size("cm")[2]
		key.length = kl/sz
	}
	if (key.pos %in% c(1, 3)) {
		ylim = c(0, 1)
		xlim = xy_from_r(range(breaks), key.length, offset)
		mar = c(0, ifelse(axes, 2.1, 1), 0, 1)
		mar[key.pos] = 2.1
	} else {
		ylim = xy_from_r(range(breaks), key.length, offset)
		xlim = c(0, 1)
		mar = c(ifelse(axes, 2.1, 1), 0, 1.2, 0)
		mar[key.pos] = ksz
	}
	par(mar = mar)

	poly = vector(mode="list", length(col))
	for (i in seq(poly))
		poly[[i]] = c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])

	tryCatch({
		plot(1, 1, type = "n", ylim = ylim, xlim = xlim, axes = FALSE,
			xlab = "", ylab = "", xaxs = "i", yaxs = "i")
		},
		error = function(x) {
			sz = cm(max(strwidth(z, "inches"))) * 1.3 + par("ps")/12 # cm
			stop(paste0("key.width too small, try key.width = lcm(", signif(sz, 3), ")"), call. = FALSE)
		}
	)

	for(i in seq_along(poly)) {
		if (key.pos %in% c(1, 3))
			polygon(poly[[i]], c(0, 0, 1, 1), col = col[i], border = NA)
		if (key.pos %in% c(2, 4))
			polygon(c(0, 0, 1, 1), poly[[i]], col = col[i], border = NA)
	}

	# box() now would draw around [0,1]:
	bx = c(breaks[1], rep(tail(breaks, 1), 2), breaks[1])
	if (key.pos %in% c(1,3))
		polygon(bx, c(0, 0, 1, 1), col = NA, border = 'black')
	if (key.pos %in% c(2,4))
		polygon(c(0, 0, 1, 1), bx, col = NA, border = 'black')

	if (add.axis) {
		opar = par(las = 1)
		axis(key.pos, at = 1:n, labels = z, cex.axis = cex.axis)
		par(opar)
	}
}

# nocov start
#' @export
identify.sfc = function(x, ..., n = min(10, length(x)), type = "n") {
	l = locator(n, type = type)
	pts = st_as_sf(as.data.frame(do.call(cbind, l)), coords = c("x", "y"), crs = st_crs(x))
	sapply(st_intersects(pts, x), function(x) if (length(x)) x[1] else NA_integer_)
}

#' @export
identify.sf = function(x, ...) {
	identify(st_geometry(x), ...)
}
# nocov end