diff --git a/NAMESPACE b/NAMESPACE index 0b5a5bc10c..6868d58b19 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -592,6 +592,7 @@ export(scale_y_sqrt) export(scale_y_time) export(sec_axis) export(set_last_plot) +export(sf_transform_xy) export(should_stop) export(stage) export(standardise_aes_names) diff --git a/NEWS.md b/NEWS.md index d750ed1592..8f6471038f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,16 @@ This is a small release focusing on fixing regressions introduced in 3.3.1. * `annotation_raster()` adds support for native rasters. For large rasters, native rasters render significantly faster than arrays (@kent37, #3388) + +* `coord_sf()` now has an argument `default_crs` that specifies the coordinate + reference system (CRS) for non-sf layers and scale/coord limits. This argument + defaults to the World Geodetic System 1984 (WGS84), which means x and y positions + are interpreted as longitude and latitude. This is a potentially breaking change + for users who use projected coordinates in non-sf layers or in limits. Setting + `default_crs = NULL` recovers the old behavior. Further, authors of extension + packages implementing `stat_sf()`-like functionality are encouraged to look at the + source code of `stat_sf()`'s `compute_group()` function to see how to provide + scale-limit hints to `coord_sf()` (@clauswilke, #3659). * Facet strips now have dedicated position-dependent theme elements (`strip.text.x.top`, `strip.text.x.bottom`, `strip.text.y.left`, diff --git a/R/annotation-map.r b/R/annotation-map.r index 73d303ec1d..8e05e2c874 100644 --- a/R/annotation-map.r +++ b/R/annotation-map.r @@ -1,34 +1,60 @@ #' @include geom-map.r NULL -#' Annotation: a maps +#' Annotation: a map #' -#' Display a fixed map on a plot. +#' Display a fixed map on a plot. This function predates the [`geom_sf()`] +#' framework and does not work with sf geometry columns as input. However, +#' it can be used in conjunction with `geom_sf()` layers and/or +#' [`coord_sf()`] (see examples). #' -#' @param map data frame representing a map. Most map objects can be -#' converted into the right format by using [fortify()] -#' @param ... other arguments used to modify aesthetics +#' @param map Data frame representing a map. See [`geom_map()`] for +#' details. +#' @param ... Other arguments used to modify visual parameters, such as +#' `colour` or `fill`. #' @export #' @examples -#' if (require("maps")) { -#' usamap <- map_data("state") +#' \dontrun{ +#' if (requireNamespace("maps", quietly = TRUE)) { +#' # location of cities in North Carolina +#' df <- data.frame( +#' name = c("Charlotte", "Raleigh", "Greensboro"), +#' lat = c(35.227, 35.772, 36.073), +#' long = c(-80.843, -78.639, -79.792) +#' ) #' -#' seal.sub <- subset(seals, long > -130 & lat < 45 & lat > 40) -#' ggplot(seal.sub, aes(x = long, y = lat)) + -#' annotation_map(usamap, fill = NA, colour = "grey50") + -#' geom_segment(aes(xend = long + delta_long, yend = lat + delta_lat)) -#' } +#' p <- ggplot(df, aes(x = long, y = lat)) + +#' annotation_map( +#' map_data("state"), +#' fill = "antiquewhite", colour = "darkgrey" +#' ) + +#' geom_point(color = "blue") + +#' geom_text( +#' aes(label = name), +#' hjust = 1.105, vjust = 1.05, color = "blue" +#' ) #' -#' if (require("maps")) { -#' seal2 <- transform(seal.sub, -#' latr = cut(lat, 2), -#' longr = cut(long, 2)) +#' # use without coord_sf() is possible but not recommended +#' p + xlim(-84, -76) + ylim(34, 37.2) #' -#' ggplot(seal2, aes(x = long, y = lat)) + -#' annotation_map(usamap, fill = NA, colour = "grey50") + -#' geom_segment(aes(xend = long + delta_long, yend = lat + delta_lat)) + -#' facet_grid(latr ~ longr, scales = "free", space = "free") -#' } +#' if (requireNamespace("sf", quietly = TRUE)) { +#' # use with coord_sf() for appropriate projection +#' p + +#' coord_sf( +#' crs = st_crs(3347), +#' xlim = c(-84, -76), +#' ylim = c(34, 37.2) +#' ) +#' +#' # you can mix annotation_map() and geom_sf() +#' nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE) +#' p + +#' geom_sf( +#' data = nc, inherit.aes = FALSE, +#' fill = NA, color = "black", size = 0.1 +#' ) + +#' coord_sf(crs = st_crs(3347)) +#' }}} annotation_map <- function(map, ...) { # Get map input into correct form if (!is.data.frame(map)) { diff --git a/R/coord-sf.R b/R/coord-sf.R index 8fc97d746f..cf99837c00 100644 --- a/R/coord-sf.R +++ b/R/coord-sf.R @@ -4,10 +4,36 @@ #' @format NULL CoordSf <- ggproto("CoordSf", CoordCartesian, - # Find the first CRS if not already supplied + # CoordSf needs to keep track of some parameters + # internally as the plot is built. These are stored + # here. + params = list(), + + # the method used to convert limits across nonlinear + # coordinate systems. + lims_method = "cross", + + get_default_crs = function(self) { + self$default_crs %||% self$params$default_crs + }, + setup_params = function(self, data) { + crs <- self$determine_crs(data) + + params <- list( + crs = crs, + default_crs = self$default_crs %||% crs + ) + self$params <- params + + params + }, + + # Helper function for setup_params(), + # finds the first CRS if not already supplied + determine_crs = function(self, data) { if (!is.null(self$crs)) { - return(list(crs = self$crs)) + return(self$crs) } for (layer_data in data) { @@ -20,10 +46,10 @@ CoordSf <- ggproto("CoordSf", CoordCartesian, if (is.na(crs)) next - return(list(crs = crs)) + return(crs) } - list(crs = NULL) + NULL }, # Transform all layers to common CRS (if provided) @@ -40,16 +66,31 @@ CoordSf <- ggproto("CoordSf", CoordCartesian, }) }, + # Allow sf layer to record the bounding boxes of elements + record_bbox = function(self, xmin, xmax, ymin, ymax) { + bbox <- self$params$bbox + bbox$xmin <- min(bbox$xmin, xmin) + bbox$xmax <- max(bbox$xmax, xmax) + bbox$ymin <- min(bbox$ymin, ymin) + bbox$ymax <- max(bbox$ymax, ymax) + self$params$bbox <- bbox + }, + transform = 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 + + # normalize geometry data, it should already be in the correct crs here data[[ geom_column(data) ]] <- sf_rescale01( data[[ geom_column(data) ]], panel_params$x_range, panel_params$y_range ) - # Assume x and y supplied directly already in common CRS + # transform and normalize regular position data data <- transform_position( - data, + sf_transform_xy(data, target_crs, source_crs), function(x) sf_rescale01_x(x, panel_params$x_range), function(x) sf_rescale01_x(x, panel_params$y_range) ) @@ -126,11 +167,51 @@ CoordSf <- ggproto("CoordSf", CoordCartesian, }, setup_panel_params = function(self, scale_x, scale_y, params = list()) { - # Bounding box of the data + # expansion factors for scale limits expansion_x <- default_expansion(scale_x, expand = self$expand) - x_range <- expand_limits_scale(scale_x, expansion_x, coord_limits = self$limits$x) expansion_y <- default_expansion(scale_y, expand = self$expand) - y_range <- expand_limits_scale(scale_y, expansion_y, coord_limits = self$limits$y) + + # get scale limits and coord limits and merge together + # coord limits take precedence over scale limits + scale_xlim <- scale_x$get_limits() + scale_ylim <- scale_y$get_limits() + coord_xlim <- self$limits$x %||% c(NA_real_, NA_real_) + coord_ylim <- self$limits$y %||% c(NA_real_, NA_real_) + + scale_xlim <- ifelse(is.na(coord_xlim), scale_xlim, coord_xlim) + scale_ylim <- ifelse(is.na(coord_ylim), scale_ylim, coord_ylim) + + # now, transform limits to common crs + # note: return value is not a true bounding box, but a + # list of x and y values whose max/mins are the bounding + # box + scales_bbox <- calc_limits_bbox( + self$lims_method, + scale_xlim, scale_ylim, + params$crs, params$default_crs + ) + + # merge coord bbox into scale limits if scale limits not explicitly set + if (is.null(self$limits$x) && is.null(self$limits$y) && + is.null(scale_x$limits) && is.null(scale_y$limits)) { + coord_bbox <- self$params$bbox + scales_xrange <- range(scales_bbox$x, coord_bbox$xmin, coord_bbox$xmax, na.rm = TRUE) + scales_yrange <- range(scales_bbox$y, coord_bbox$ymin, coord_bbox$ymax, na.rm = TRUE) + } else if (any(!is.finite(scales_bbox$x) | !is.finite(scales_bbox$y))) { + if (self$lims_method != "geometry_bbox") { + warn("Projection of x or y limits failed in `coord_sf()`.\nConsider setting `lims_method = \"geometry_bbox\"` or `default_crs = NULL`.") + } + coord_bbox <- self$params$bbox + scales_xrange <- c(coord_bbox$xmin, coord_bbox$xmax) %||% c(0, 0) + scales_yrange <- c(coord_bbox$ymin, coord_bbox$ymax) %||% c(0, 0) + } else { + scales_xrange <- range(scales_bbox$x, na.rm = TRUE) + scales_yrange <- range(scales_bbox$y, na.rm = TRUE) + } + + # apply coordinate expansion + x_range <- expand_limits_continuous(scales_xrange, expansion_x) + y_range <- expand_limits_continuous(scales_yrange, expansion_y) bbox <- c( x_range[1], y_range[1], x_range[2], y_range[2] @@ -160,25 +241,39 @@ CoordSf <- ggproto("CoordSf", CoordCartesian, y_range = y_range, graticule = graticule, crs = params$crs, + default_crs = params$default_crs, label_axes = self$label_axes, label_graticule = self$label_graticule ) }, - backtransform_range = function(panel_params) { - # this does not actually return backtransformed ranges in the general case, needs fixing - warn("range backtransformation not implemented in this coord; results may be wrong.") - list(x = panel_params$x_range, y = panel_params$y_range) + backtransform_range = function(self, panel_params) { + target_crs <- panel_params$default_crs + source_crs <- panel_params$crs + + x <- panel_params$x_range + y <- panel_params$y_range + data <- list(x = c(x, x), y = c(y, rev(y))) + data <- sf_transform_xy(data, target_crs, source_crs) + list(x = range(data$x), y = range(data$y)) }, range = function(panel_params) { list(x = panel_params$x_range, y = panel_params$y_range) }, - # CoordSf enforces a fixed aspect ratio -> axes cannot be changed freely under faceting is_free = function() FALSE, + # for regular geoms (such as geom_path, geom_polygon, etc.), CoordSf is non-linear + is_linear = function() FALSE, + + distance = function(self, x, y, panel_params) { + d <- self$backtransform_range(panel_params) + max_dist <- dist_euclidean(d$x, d$y) + dist_euclidean(x, y) / max_dist + }, + aspect = function(self, panel_params) { if (isTRUE(sf::st_is_longlat(panel_params$crs))) { # Contributed by @edzer @@ -372,6 +467,58 @@ CoordSf <- ggproto("CoordSf", CoordCartesian, } ) +#' Transform spatial position data +#' +#' Helper function that can transform spatial position data (pairs of x, y +#' values) among coordinate systems. +#' +#' @param data Data frame or list containing numerical columns `x` and `y`. +#' @param target_crs,source_crs Target and source coordinate reference systems. +#' If `NULL` or `NA`, the data is not transformed. +#' @return A copy of the input data with `x` and `y` replaced by transformed values. +#' @examples +#' if (requireNamespace("sf", quietly = TRUE)) { +#' # location of cities in NC by long (x) and lat (y) +#' data <- data.frame( +#' city = c("Charlotte", "Raleigh", "Greensboro"), +#' x = c(-80.843, -78.639, -79.792), +#' y = c(35.227, 35.772, 36.073) +#' ) +#' +#' # transform to projected coordinates +#' data_proj <- sf_transform_xy(data, 3347, 4326) +#' data_proj +#' +#' # transform back +#' sf_transform_xy(data_proj, 4326, 3347) +#' } +#' @keywords internal +#' @export +sf_transform_xy <- function(data, target_crs, source_crs) { + if (identical(target_crs, source_crs) || + is.null(target_crs) || is.null(source_crs) || is.null(data) || + is.na(target_crs) || is.na(source_crs) || + !all(c("x", "y") %in% names(data))) { + return(data) + } + + # by turning the data into a geometry list column of individual points, + # we can make sure that the output length equals the input length, even + # if the transformation fails in some cases + sf_data <- sf::st_sfc( + mapply(function(x, y) sf::st_point(as.numeric(c(x, y))), data$x, data$y, SIMPLIFY = FALSE), + crs = source_crs + ) + sf_data_trans <- sf::st_transform(sf_data, target_crs) + data$x <- vapply(sf_data_trans, function(x) x[1], numeric(1)) + data$y <- vapply(sf_data_trans, function(x) x[2], numeric(1)) + + data +} + +## helper functions to normalize geometry and position data + +# normalize geometry data (variable x is geometry column) sf_rescale01 <- function(x, x_range, y_range) { if (is.null(x)) { return(x) @@ -379,14 +526,98 @@ sf_rescale01 <- function(x, x_range, y_range) { sf::st_normalize(x, c(x_range[1], y_range[1], x_range[2], y_range[2])) } + +# normalize position data (variable x is x or y position) sf_rescale01_x <- function(x, range) { (x - range[1]) / diff(range) } +# different limits methods +calc_limits_bbox <- function(method, xlim, ylim, crs, default_crs) { + if (any(!is.finite(c(xlim, ylim))) && method != "geometry_bbox") { + abort("Scale limits cannot be mapped onto spatial coordinates in `coord_sf()`.\nConsider setting `lims_method = \"geometry_bbox\"` or `default_crs = NULL`.") + } + + bbox <- switch( + method, + # For method "box", we take the limits and turn them into a + # box. We subdivide the box edges into multiple segments to + # better cover the respective area under non-linear transformation + box = list( + x = c( + rep(xlim[1], 20), seq(xlim[1], xlim[2], length.out = 20), + rep(xlim[2], 20), seq(xlim[2], xlim[1], length.out = 20) + ), + y = c( + seq(ylim[1], ylim[2], length.out = 20), rep(ylim[2], 20), + seq(ylim[2], ylim[1], length.out = 20), rep(ylim[1], 20) + ) + ), + # For method "geometry_bbox" we ignore all limits info provided here + geometry_bbox = list( + x = c(NA_real_, NA_real_), + y = c(NA_real_, NA_real_) + ), + # For method "orthogonal" we simply return what we are given + orthogonal = list( + x = xlim, + y = ylim + ), + # For method "cross" we take the mid-point along each side of + # the scale range for better behavior when box is nonlinear or + # rotated in projected space + # + # Method "cross" is also the default + cross =, + list( + x = c(rep(mean(xlim), 20), seq(xlim[1], xlim[2], length.out = 20)), + y = c(seq(ylim[1], ylim[2], length.out = 20), rep(mean(ylim), 20)) + ) + ) + sf_transform_xy(bbox, crs, default_crs) +} -#' @param crs Use this to select a specific coordinate reference system (CRS). -#' If not specified, will use the CRS defined in the first layer. -#' @param datum CRS that provides datum to use when generating graticules +#' @param crs The coordinate reference system (CRS) into which all data should +#' be projected before plotting. If not specified, will use the CRS defined +#' in the first sf layer of the plot. +#' @param default_crs The default CRS to be used for non-sf layers (which +#' don't carry any CRS information) and scale limits. If not specified, this +#' defaults to the World Geodetic System 1984 (WGS84), which means x and y +#' positions are interpreted as longitude and latitude, respectively. If +#' set to `NULL`, uses the setting for `crs`, which means that then all +#' non-sf layers and scale limits are assumed to be specified in projected +#' coordinates. +#' @param xlim,ylim Limits for the x and y axes. These limits are specified +#' in the units of the default CRS. To specify limits in projected coordinates, +#' set `default_crs = NULL`. How limit specifications translate into the exact +#' region shown on the plot can be confusing when non-linear or rotated coordinate +#' systems are used. First, different methods can be preferable under different +#' conditions. See parameter `lims_method` for details. Second, specifying limits +#' along only one direction can affect the automatically generated limits along the +#' other direction. Therefore, it is best to always specify limits for both x and y. +#' Third, specifying limits via position scales or `xlim()`/`ylim()` is strongly +#' discouraged, as it can result in data points being dropped from the plot even +#' though they would be visible in the final plot region. Finally, specifying limits +#' that cross the international date boundary is not possible with WGS84 as the default +#' crs. All these issues can be avoided by working in projected coordinates, +#' via `default_crs = NULL`, but at the cost of having to provide less intuitive +#' numeric values for the limit parameters. +#' @param lims_method Method specifying how scale limits are converted into +#' limits on the plot region. For a very non-linear CRS (e.g., a perspective centered +#' around the North pole), the available methods yield widely differing results, and +#' you may want to try various options. Methods currently implemented include `"cross"` +#' (the default), `"box"`, `"orthogonal"`, and `"geometry_bbox"`. For method `"cross"`, +#' limits along one direction (e.g., longitude) are applied at the midpoint of the +#' other direction (e.g., latitude). This method avoids excessively large limits for +#' rotated coordinate systems but means that sometimes limits need to be expanded a +#' little further if extreme data points are to be included in the final plot region. +#' By contrast, for method `"box"`, a box is generated out of the limits along both directions, +#' and then limits in projected coordinates are chosen such that the entire box is +#' visible. This method can yield plot regions that are too large. Finally, method +#' `"orthogonal"` applies limits separately along each axis, and method +#' `"geometry_bbox"` ignores all limit information except the bounding boxes of any +#' objects in the `geometry` aesthetic. +#' @param datum CRS that provides datum to use when generating graticules. #' @param label_axes Character vector or named list of character values #' specifying which graticule lines (meridians or parallels) should be labeled on #' which side of the plot. Meridians are indicated by `"E"` (for East) and @@ -411,15 +642,16 @@ sf_rescale01_x <- function(x, range) { #' generally yield better results and should be used instead. #' #' This parameter can be used alone or in combination with `label_axes`. -#' @param ndiscr number of segments to use for discretising graticule lines; -#' try increasing this when graticules look unexpected +#' @param ndiscr Number of segments to use for discretising graticule lines; +#' try increasing this number when graticules look incorrect. #' @inheritParams coord_cartesian #' @export #' @rdname ggsf coord_sf <- function(xlim = NULL, ylim = NULL, expand = TRUE, - crs = NULL, datum = sf::st_crs(4326), + crs = NULL, default_crs = sf::st_crs(4326), + datum = sf::st_crs(4326), label_graticule = waiver(), - label_axes = waiver(), + label_axes = waiver(), lims_method = c("cross", "box", "orthogonal", "geometry_bbox"), ndiscr = 100, default = FALSE, clip = "on") { if (is.waive(label_graticule) && is.waive(label_axes)) { @@ -447,10 +679,19 @@ coord_sf <- function(xlim = NULL, ylim = NULL, expand = TRUE, label_graticule <- "" } + # switch limit method to "orthogonal" if not specified and default_crs indicates projected coords + if (is.null(default_crs) && is_missing(lims_method)) { + lims_method <- "orthogonal" + } else { + lims_method <- match.arg(lims_method) + } + ggproto(NULL, CoordSf, limits = list(x = xlim, y = ylim), + lims_method = lims_method, datum = datum, crs = crs, + default_crs = default_crs, label_axes = label_axes, label_graticule = label_graticule, ndiscr = ndiscr, diff --git a/R/geom-sf.R b/R/geom-sf.R index 56e156f1aa..ac8988fc09 100644 --- a/R/geom-sf.R +++ b/R/geom-sf.R @@ -27,6 +27,29 @@ #' either specify it using the `CRS` param, or `coord_sf()` will #' take it from the first layer that defines a CRS. #' +#' @section Combining sf layers and regular geoms: +#' Most regular geoms, such as [geom_point()], [geom_path()], +#' [geom_text()], [geom_polygon()] etc. will work fine with `coord_sf()`. However +#' when using these geoms, two problems arise. First, what CRS should be used +#' for the x and y coordinates used by these non-sf geoms? The CRS applied to +#' non-sf geoms is set by the `default_crs` parameter, and it defaults to +#' the World Geodetic System 1984 (WGS84). This means that x and y +#' positions are interpreted as longitude and latitude, respectively. You can +#' also specify any other valid CRS as the default CRS for non-sf geoms. Moreover, +#' if you set `default_crs = NULL`, then positions for non-sf geoms are +#' interpreted as projected coordinates. This setting allows you complete control +#' over where exactly items are placed on the plot canvas. +#' +#' The second problem that arises for non-sf geoms is how straight lines +#' should be interpreted in projected space. The approach `coord_sf()` takes is +#' to break straight lines into small pieces (i.e., segmentize them) and +#' then transform the pieces into projected coordinates. For the default setting +#' where x and y are interpreted as longitude and latitude, this approach means +#' that horizontal lines follow the parallels and vertical lines follow the +#' meridians. If you need a different approach to handling straight lines, then +#' you should manually segmentize and project coordinates and generate the plot +#' in projected coordinates. +#' #' @param show.legend logical. Should this layer be included in the legends? #' `NA`, the default, includes if any aesthetics are mapped. #' `FALSE` never includes, and `TRUE` always includes. diff --git a/R/stat-sf-coordinates.R b/R/stat-sf-coordinates.R index f1616ce24a..86daf026b0 100644 --- a/R/stat-sf-coordinates.R +++ b/R/stat-sf-coordinates.R @@ -84,12 +84,36 @@ stat_sf_coordinates <- function(mapping = aes(), data = NULL, geom = "point", #' @export StatSfCoordinates <- ggproto( "StatSfCoordinates", Stat, - compute_group = function(data, scales, fun.geometry = NULL) { + + compute_layer = function(self, data, params, layout) { + # add coord to the params, so it can be forwarded to compute_group() + params$coord <- layout$coord + ggproto_parent(Stat, self)$compute_layer(data, params, layout) + }, + + compute_group = function(self, data, scales, coord, fun.geometry = NULL) { if (is.null(fun.geometry)) { fun.geometry <- function(x) sf::st_point_on_surface(sf::st_zm(x)) } points_sfc <- fun.geometry(data$geometry) + + if (inherits(coord, "CoordSf")) { + # register bounding box if the coord derives from CoordSf + bbox <- sf::st_bbox(points_sfc) + + coord$record_bbox( + xmin = bbox[["xmin"]], xmax = bbox[["xmax"]], + ymin = bbox[["ymin"]], ymax = bbox[["ymax"]] + ) + + # transform to the coord's default crs if possible + default_crs <- coord$get_default_crs() + if (!(is.null(default_crs) || is.na(default_crs) || + is.na(sf::st_crs(points_sfc)))) { + points_sfc <- sf::st_transform(points_sfc, default_crs) + } + } coordinates <- sf::st_coordinates(points_sfc) data$x <- coordinates[, "X"] data$y <- coordinates[, "Y"] diff --git a/R/stat-sf.R b/R/stat-sf.R index 3e2a840d78..9e8e139f77 100644 --- a/R/stat-sf.R +++ b/R/stat-sf.R @@ -3,12 +3,52 @@ #' @usage NULL #' @format NULL StatSf <- ggproto("StatSf", Stat, - compute_group = function(data, scales) { - bbox <- sf::st_bbox(data[[ geom_column(data) ]]) - data$xmin <- bbox[["xmin"]] - data$xmax <- bbox[["xmax"]] - data$ymin <- bbox[["ymin"]] - data$ymax <- bbox[["ymax"]] + compute_layer = function(self, data, params, layout) { + # add coord to the params, so it can be forwarded to compute_group() + params$coord <- layout$coord + ggproto_parent(Stat, self)$compute_layer(data, params, layout) + }, + + compute_group = function(data, scales, coord) { + geometry_data <- data[[ geom_column(data) ]] + geometry_crs <- sf::st_crs(geometry_data) + + bbox <- sf::st_bbox(geometry_data) + + if (inherits(coord, "CoordSf")) { + # if the coord derives from CoordSf, then it + # needs to know about bounding boxes of geometry data + coord$record_bbox( + xmin = bbox[["xmin"]], xmax = bbox[["xmax"]], + ymin = bbox[["ymin"]], ymax = bbox[["ymax"]] + ) + + # to represent the location of the geometry in default coordinates, + # we take the mid-point along each side of the bounding box and + # backtransform + bbox_trans <- sf_transform_xy( + list( + x = c(rep(0.5*(bbox[["xmin"]] + bbox[["xmax"]]), 2), bbox[["xmin"]], bbox[["xmax"]]), + y = c(bbox[["ymin"]], bbox[["ymax"]], rep(0.5*(bbox[["ymin"]] + bbox[["ymax"]]), 2)) + ), + coord$get_default_crs(), + geometry_crs + ) + + # record as xmin, xmax, ymin, ymax so regular scales + # have some indication of where shapes lie + data$xmin <- min(bbox_trans$x) + data$xmax <- max(bbox_trans$x) + data$ymin <- min(bbox_trans$y) + data$ymax <- max(bbox_trans$y) + } else { + # for all other coords, we record the full extent of the + # geometry object + data$xmin <- bbox[["xmin"]] + data$xmax <- bbox[["xmax"]] + data$ymin <- bbox[["ymin"]] + data$ymax <- bbox[["ymax"]] + } data }, diff --git a/man/annotation_map.Rd b/man/annotation_map.Rd index b093f0be03..468bd4c049 100644 --- a/man/annotation_map.Rd +++ b/man/annotation_map.Rd @@ -2,37 +2,63 @@ % Please edit documentation in R/annotation-map.r \name{annotation_map} \alias{annotation_map} -\title{Annotation: a maps} +\title{Annotation: a map} \usage{ annotation_map(map, ...) } \arguments{ -\item{map}{data frame representing a map. Most map objects can be -converted into the right format by using \code{\link[=fortify]{fortify()}}} +\item{map}{Data frame representing a map. See \code{\link[=geom_map]{geom_map()}} for +details.} -\item{...}{other arguments used to modify aesthetics} +\item{...}{Other arguments used to modify visual parameters, such as +\code{colour} or \code{fill}.} } \description{ -Display a fixed map on a plot. +Display a fixed map on a plot. This function predates the \code{\link[=geom_sf]{geom_sf()}} +framework and does not work with sf geometry columns as input. However, +it can be used in conjunction with \code{geom_sf()} layers and/or +\code{\link[=coord_sf]{coord_sf()}} (see examples). } \examples{ -if (require("maps")) { -usamap <- map_data("state") +\dontrun{ +if (requireNamespace("maps", quietly = TRUE)) { +# location of cities in North Carolina +df <- data.frame( + name = c("Charlotte", "Raleigh", "Greensboro"), + lat = c(35.227, 35.772, 36.073), + long = c(-80.843, -78.639, -79.792) +) -seal.sub <- subset(seals, long > -130 & lat < 45 & lat > 40) -ggplot(seal.sub, aes(x = long, y = lat)) + - annotation_map(usamap, fill = NA, colour = "grey50") + - geom_segment(aes(xend = long + delta_long, yend = lat + delta_lat)) -} +p <- ggplot(df, aes(x = long, y = lat)) + + annotation_map( + map_data("state"), + fill = "antiquewhite", colour = "darkgrey" + ) + + geom_point(color = "blue") + + geom_text( + aes(label = name), + hjust = 1.105, vjust = 1.05, color = "blue" + ) -if (require("maps")) { -seal2 <- transform(seal.sub, - latr = cut(lat, 2), - longr = cut(long, 2)) +# use without coord_sf() is possible but not recommended +p + xlim(-84, -76) + ylim(34, 37.2) -ggplot(seal2, aes(x = long, y = lat)) + - annotation_map(usamap, fill = NA, colour = "grey50") + - geom_segment(aes(xend = long + delta_long, yend = lat + delta_lat)) + - facet_grid(latr ~ longr, scales = "free", space = "free") -} +if (requireNamespace("sf", quietly = TRUE)) { +# use with coord_sf() for appropriate projection +p + + coord_sf( + crs = st_crs(3347), + xlim = c(-84, -76), + ylim = c(34, 37.2) + ) + +# you can mix annotation_map() and geom_sf() +nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE) +p + + geom_sf( + data = nc, inherit.aes = FALSE, + fill = NA, color = "black", size = 0.1 + ) + + coord_sf(crs = st_crs(3347)) +}}} } diff --git a/man/ggsf.Rd b/man/ggsf.Rd index cdfdb88b97..93c15928c6 100644 --- a/man/ggsf.Rd +++ b/man/ggsf.Rd @@ -18,9 +18,11 @@ coord_sf( ylim = NULL, expand = TRUE, crs = NULL, + default_crs = sf::st_crs(4326), datum = sf::st_crs(4326), label_graticule = waiver(), label_axes = waiver(), + lims_method = c("cross", "box", "orthogonal", "geometry_bbox"), ndiscr = 100, default = FALSE, clip = "on" @@ -83,18 +85,39 @@ stat_sf( ) } \arguments{ -\item{xlim}{Limits for the x and y axes.} - -\item{ylim}{Limits for the x and y axes.} +\item{xlim, ylim}{Limits for the x and y axes. These limits are specified +in the units of the default CRS. To specify limits in projected coordinates, +set \code{default_crs = NULL}. How limit specifications translate into the exact +region shown on the plot can be confusing when non-linear or rotated coordinate +systems are used. First, different methods can be preferable under different +conditions. See parameter \code{lims_method} for details. Second, specifying limits +along only one direction can affect the automatically generated limits along the +other direction. Therefore, it is best to always specify limits for both x and y. +Third, specifying limits via position scales or \code{xlim()}/\code{ylim()} is strongly +discouraged, as it can result in data points being dropped from the plot even +though they would be visible in the final plot region. Finally, specifying limits +that cross the international date boundary is not possible with WGS84 as the default +crs. All these issues can be avoided by working in projected coordinates, +via \code{default_crs = NULL}, but at the cost of having to provide less intuitive +numeric values for the limit parameters.} \item{expand}{If \code{TRUE}, the default, adds a small expansion factor to the limits to ensure that data and axes don't overlap. If \code{FALSE}, limits are taken exactly from the data or \code{xlim}/\code{ylim}.} -\item{crs}{Use this to select a specific coordinate reference system (CRS). -If not specified, will use the CRS defined in the first layer.} +\item{crs}{The coordinate reference system (CRS) into which all data should +be projected before plotting. If not specified, will use the CRS defined +in the first sf layer of the plot.} + +\item{default_crs}{The default CRS to be used for non-sf layers (which +don't carry any CRS information) and scale limits. If not specified, this +defaults to the World Geodetic System 1984 (WGS84), which means x and y +positions are interpreted as longitude and latitude, respectively. If +set to \code{NULL}, uses the setting for \code{crs}, which means that then all +non-sf layers and scale limits are assumed to be specified in projected +coordinates.} -\item{datum}{CRS that provides datum to use when generating graticules} +\item{datum}{CRS that provides datum to use when generating graticules.} \item{label_graticule}{Character vector indicating which graticule lines should be labeled where. Meridians run north-south, and the letters \code{"N"} and \code{"S"} indicate that @@ -122,8 +145,24 @@ specified with \code{list(bottom = "E", left = "N")}. This parameter can be used alone or in combination with \code{label_graticule}.} -\item{ndiscr}{number of segments to use for discretising graticule lines; -try increasing this when graticules look unexpected} +\item{lims_method}{Method specifying how scale limits are converted into +limits on the plot region. For a very non-linear CRS (e.g., a perspective centered +around the North pole), the available methods yield widely differing results, and +you may want to try various options. Methods currently implemented include \code{"cross"} +(the default), \code{"box"}, \code{"orthogonal"}, and \code{"geometry_bbox"}. For method \code{"cross"}, +limits along one direction (e.g., longitude) are applied at the midpoint of the +other direction (e.g., latitude). This method avoids excessively large limits for +rotated coordinate systems but means that sometimes limits need to be expanded a +little further if extreme data points are to be included in the final plot region. +By contrast, for method \code{"box"}, a box is generated out of the limits along both directions, +and then limits in projected coordinates are chosen such that the entire box is +visible. This method can yield plot regions that are too large. Finally, method +\code{"orthogonal"} applies limits separately along each axis, and method +\code{"geometry_bbox"} ignores all limit information except the bounding boxes of any +objects in the \code{geometry} aesthetic.} + +\item{ndiscr}{Number of segments to use for discretising graticule lines; +try increasing this number when graticules look incorrect.} \item{default}{Is this the default coordinate system? If \code{FALSE} (the default), then replacing this coordinate system with another one creates a message alerting @@ -248,6 +287,31 @@ either specify it using the \code{CRS} param, or \code{coord_sf()} will take it from the first layer that defines a CRS. } +\section{Combining sf layers and regular geoms}{ + +Most regular geoms, such as \code{\link[=geom_point]{geom_point()}}, \code{\link[=geom_path]{geom_path()}}, +\code{\link[=geom_text]{geom_text()}}, \code{\link[=geom_polygon]{geom_polygon()}} etc. will work fine with \code{coord_sf()}. However +when using these geoms, two problems arise. First, what CRS should be used +for the x and y coordinates used by these non-sf geoms? The CRS applied to +non-sf geoms is set by the \code{default_crs} parameter, and it defaults to +the World Geodetic System 1984 (WGS84). This means that x and y +positions are interpreted as longitude and latitude, respectively. You can +also specify any other valid CRS as the default CRS for non-sf geoms. Moreover, +if you set \code{default_crs = NULL}, then positions for non-sf geoms are +interpreted as projected coordinates. This setting allows you complete control +over where exactly items are placed on the plot canvas. + +The second problem that arises for non-sf geoms is how straight lines +should be interpreted in projected space. The approach \code{coord_sf()} takes is +to break straight lines into small pieces (i.e., segmentize them) and +then transform the pieces into projected coordinates. For the default setting +where x and y are interpreted as longitude and latitude, this approach means +that horizontal lines follow the parallels and vertical lines follow the +meridians. If you need a different approach to handling straight lines, then +you should manually segmentize and project coordinates and generate the plot +in projected coordinates. +} + \examples{ if (requireNamespace("sf", quietly = TRUE)) { nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE) diff --git a/man/sf_transform_xy.Rd b/man/sf_transform_xy.Rd new file mode 100644 index 0000000000..3db14eb8cb --- /dev/null +++ b/man/sf_transform_xy.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/coord-sf.R +\name{sf_transform_xy} +\alias{sf_transform_xy} +\title{Transform spatial position data} +\usage{ +sf_transform_xy(data, target_crs, source_crs) +} +\arguments{ +\item{data}{Data frame or list containing numerical columns \code{x} and \code{y}.} + +\item{target_crs, source_crs}{Target and source coordinate reference systems. +If \code{NULL} or \code{NA}, the data is not transformed.} +} +\value{ +A copy of the input data with \code{x} and \code{y} replaced by transformed values. +} +\description{ +Helper function that can transform spatial position data (pairs of x, y +values) among coordinate systems. +} +\examples{ +if (requireNamespace("sf", quietly = TRUE)) { +# location of cities in NC by long (x) and lat (y) +data <- data.frame( + city = c("Charlotte", "Raleigh", "Greensboro"), + x = c(-80.843, -78.639, -79.792), + y = c(35.227, 35.772, 36.073) +) + +# transform to projected coordinates +data_proj <- sf_transform_xy(data, 3347, 4326) +data_proj + +# transform back +sf_transform_xy(data_proj, 4326, 3347) +} +} +\keyword{internal} diff --git a/tests/figs/coord-sf/default-crs-turned-off.svg b/tests/figs/coord-sf/default-crs-turned-off.svg new file mode 100644 index 0000000000..a0923bdb7b --- /dev/null +++ b/tests/figs/coord-sf/default-crs-turned-off.svg @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +35 +° +N +36 +° +N +37 +° +N +38 +° +N +39 +° +N +40 +° +N + + + + + + + + + + + +81 +° +W +80 +° +W +79 +° +W +78 +° +W +77 +° +W +x +y +default crs turned off + diff --git a/tests/figs/coord-sf/limits-specified-in-long-lat.svg b/tests/figs/coord-sf/limits-specified-in-long-lat.svg new file mode 100644 index 0000000000..e128364288 --- /dev/null +++ b/tests/figs/coord-sf/limits-specified-in-long-lat.svg @@ -0,0 +1,85 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +37 +° +N +38 +° +N +39 +° +N +40 +° +N +41 +° +N + + + + + + + + + + +81 +° +W +80 +° +W +79 +° +W +78 +° +W +77 +° +W +x +y +limits specified in long-lat + diff --git a/tests/figs/coord-sf/limits-specified-in-projected-coords.svg b/tests/figs/coord-sf/limits-specified-in-projected-coords.svg new file mode 100644 index 0000000000..9a46fdabc0 --- /dev/null +++ b/tests/figs/coord-sf/limits-specified-in-projected-coords.svg @@ -0,0 +1,85 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +37 +° +N +38 +° +N +39 +° +N +40 +° +N +41 +° +N + + + + + + + + + + +81 +° +W +80 +° +W +79 +° +W +78 +° +W +77 +° +W +x +y +limits specified in projected coords + diff --git a/tests/figs/coord-sf/non-sf-geoms-use-long-lat.svg b/tests/figs/coord-sf/non-sf-geoms-use-long-lat.svg new file mode 100644 index 0000000000..ce1f75cb76 --- /dev/null +++ b/tests/figs/coord-sf/non-sf-geoms-use-long-lat.svg @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +35 +° +N +36 +° +N +37 +° +N +38 +° +N +39 +° +N +40 +° +N + + + + + + + + + + + +81 +° +W +80 +° +W +79 +° +W +78 +° +W +77 +° +W +x +y +non-sf geoms use long-lat + diff --git a/tests/testthat/test-coord_sf.R b/tests/testthat/test-coord_sf.R index cf4d325b25..b2eac12d2f 100644 --- a/tests/testthat/test-coord_sf.R +++ b/tests/testthat/test-coord_sf.R @@ -204,3 +204,81 @@ test_that("Inf is squished to range", { expect_equal(d[[2]]$x, 0) expect_equal(d[[2]]$y, 1) }) + +test_that("default crs works", { + skip_if_not_installed("sf") + + polygon <- sf::st_sfc( + sf::st_polygon(list(matrix(c(-80, -76, -76, -80, -80, 35, 35, 40, 40, 35), ncol = 2))), + crs = 4326 # basic long-lat crs + ) + polygon <- sf::st_transform(polygon, crs = 3347) + + points <- data_frame( + x = c(-80, -80, -76, -76), + y = c(35, 40, 35, 40) + ) + + p <- ggplot(polygon) + geom_sf(fill = NA) + + # projected sf objects can be mixed with regular geoms using non-projected data + expect_doppelganger( + "non-sf geoms use long-lat", + p + geom_point(data = points, aes(x, y)) + ) + + # default crs can be turned off + points_trans <- sf_transform_xy(points, 3347, 4326) + expect_doppelganger( + "default crs turned off", + p + geom_point(data = points_trans, aes(x, y)) + + coord_sf(default_crs = NULL) + ) + + # by default, coord limits are specified in long-lat + expect_doppelganger( + "limits specified in long-lat", + p + geom_point(data = points, aes(x, y)) + + coord_sf(xlim = c(-80.5, -76), ylim = c(36, 41)) + ) + + # when default crs is off, limits are specified in projected coords + lims <- sf_transform_xy( + list(x = c(-80.5, -76, -78.25, -78.25), y = c(38.5, 38.5, 36, 41)), + 3347, 4326 + ) + expect_doppelganger( + "limits specified in projected coords", + p + geom_point(data = points_trans, aes(x, y)) + + coord_sf(xlim = lims$x[1:2], ylim = lims$y[3:4], default_crs = NULL) + ) +}) + +test_that("sf_transform_xy() works", { + skip_if_not_installed("sf") + + data <- list( + city = c("Charlotte", "Raleigh", "Greensboro"), + x = c(-80.843, -78.639, -79.792), + y = c(35.227, 35.772, 36.073) + ) + + # no transformation if one crs is missing + out <- sf_transform_xy(data, NULL, 4326) + expect_identical(data, out) + out <- sf_transform_xy(data, 4326, NULL) + expect_identical(data, out) + + # transform to projected coordinates + out <- sf_transform_xy(data, 3347, 4326) + expect_identical(data$city, out$city) # columns other than x, y are not changed + expect_true(all(abs(out$x - c(7275499, 7474260, 7357835)) < 10)) + expect_true(all(abs(out$y - c(-60169, 44384, 57438)) < 10)) + + # transform back + out2 <- sf_transform_xy(out, 4326, 3347) + expect_identical(data$city, out2$city) + expect_true(all(abs(out2$x - data$x) < .01)) + expect_true(all(abs(out2$y - data$y) < .01)) + +})