Implement default crs for non-sf objects in coord_sf().

DESCRIPTION

@@ -253,6 +253,6 @@ Collate:
     'zxx.r'
     'zzz.r'
 VignetteBuilder: knitr
-RoxygenNote: 7.0.1
+RoxygenNote: 7.0.2
 Roxygen: list(markdown = TRUE)
 Encoding: UTF-8

R/coord-sf.R

@@ -4,10 +4,32 @@
 #' @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(),
+
+  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 +42,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 +62,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 +163,34 @@ 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 transform to common crs
+    scale_xlim <- scale_x$get_limits()
+    scale_ylim <- scale_y$get_limits()
+
+    scales_bbox <- sf_transform_xy(
+      list(x = c(scale_xlim, scale_xlim), y = c(scale_ylim, rev(scale_ylim))),
+      params$crs, params$default_crs
+    )
+
+    # merge scale limits and coord limits
+    coord_bbox <- self$params$bbox
+    scales_xrange <- c(min(scales_bbox$x, coord_bbox$xmin), max(scales_bbox$x, coord_bbox$xmax))
+    scales_yrange <- c(min(scales_bbox$y, coord_bbox$ymin), max(scales_bbox$y, coord_bbox$ymax))
+
+    # calculate final coord limits by putting everything together and applying expansion
+    coord_limits_x <- self$limits$x %||% c(NA_real_, NA_real_)
+    coord_limits_y <- self$limits$y %||% c(NA_real_, NA_real_)
+
+    x_range <- expand_limits_continuous(
+      scales_xrange, expansion_x, coord_limits = coord_limits_x
+    )
+    y_range <- expand_limits_continuous(
+      scales_yrange, expansion_y, coord_limits = coord_limits_y
+    )
     bbox <- c(
       x_range[1], y_range[1],
       x_range[2], y_range[2]
@@ -160,28 +220,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
-    warning(
-      "range backtransformation not implemented in this coord; results may be wrong.",
-      call. = FALSE
-    )
-    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
@@ -375,20 +446,63 @@ CoordSf <- ggproto("CoordSf", CoordCartesian,
   }
 )
 
+## helper functions to transform and normalize geometry and position data
+# transform position data (columns x and y in a data frame)
+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)
+  }
+
+  # we need to exclude any non-finite values from st_transform
+  # we replace them with 0 and afterwards with NA
+  finite_x <- is.finite(data$x)
+  finite_y <- is.finite(data$y)
+  data$x[!finite_x] <- 0
+  data$y[!finite_y] <- 0
+
+  sf_data <- sf::st_sfc(
+    sf::st_multipoint(cbind(data$x, data$y)),
+    crs = source_crs
+  )
+  sf_data_trans <- sf::st_transform(sf_data, target_crs)[[1]]
+  data$x <- sf_data_trans[, 1]
+  data$y <- sf_data_trans[, 2]
+
+  data$x[!(finite_x & finite_y)] <- NA
+  data$y[!(finite_x & finite_y)] <- NA
+  data
+}
+
+# normalize geometry data (variable x is geometry column)
 sf_rescale01 <- function(x, x_range, y_range) {
   if (is.null(x)) {
     return(x)
   }
 
   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)
 }
 
 
+
 #' @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 default_crs The default CRS to be used for non-sf layers (which
+#'   don't carry any CRS information). 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. The default CRS
+#'   is also the reference system used to set limits via position scales. If
+#'   set to `NULL`, uses the setting for `crs`.
+#' @param xlim,ylim Limits for the x and y axes. These limits are specified
+#'   in the units of the CRS set via the `crs` argument or, if `crs` is not
+#'   specified, the CRS of the first layer that has a CRS.
 #' @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
@@ -417,7 +531,8 @@ sf_rescale01_x <- function(x, range) {
 #' @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(),
                      ndiscr = 100, default = FALSE, clip = "on") {
@@ -457,6 +572,7 @@ coord_sf <- function(xlim = NULL, ylim = NULL, expand = TRUE,
     limits = list(x = xlim, y = ylim),
     datum = datum,
     crs = crs,
+    default_crs = default_crs,
     label_axes = label_axes,
     label_graticule = label_graticule,
     ndiscr = ndiscr,

R/stat-sf-coordinates.R

@@ -84,12 +84,28 @@ 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)
+
+    # transform to the coord's default crs if possible
+    if (inherits(coord, "CoordSf")) {
+      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"]

R/stat-sf.R

@@ -3,12 +3,49 @@
 #' @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"]]
+      )
+
+      # register geometric center of each bbox, to give regular scales
+      # some indication of where shapes lie
+      bbox_trans <- sf_transform_xy(
+        list(
+          x = 0.5*(bbox[["xmin"]] + bbox[["xmax"]]),
+          y = 0.5*(bbox[["ymin"]] + bbox[["ymax"]])
+        ),
+        coord$get_default_crs(),
+        geometry_crs
+      )
+
+      data$xmin <- bbox_trans$x
+      data$xmax <- bbox_trans$x
+      data$ymin <- bbox_trans$y
+      data$ymax <- 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
   },