1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
|
#' Map projections
#'
#' @description
#' `r lifecycle::badge("superseded")`
#'
#' `coord_map()` projects a portion of the earth, which is approximately
#' spherical, onto a flat 2D plane using any projection defined by the
#' `mapproj` package. Map projections do not, in general, preserve straight
#' lines, so this requires considerable computation. `coord_quickmap()` is a
#' quick approximation that does preserve straight lines. It works best for
#' smaller areas closer to the equator.
#'
#' Both `coord_map()` and `coord_quickmap()`
#' are superseded by [`coord_sf()`], and should no longer be used in new
#' code. All regular (non-sf) geoms can be used with `coord_sf()` by
#' setting the default coordinate system via the `default_crs` argument.
#' See also the examples for [`annotation_map()`] and [`geom_map()`].
#'
#' @details
#'
#' Map projections must account for the fact that the actual length
#' (in km) of one degree of longitude varies between the equator and the pole.
#' Near the equator, the ratio between the lengths of one degree of latitude and
#' one degree of longitude is approximately 1. Near the pole, it tends
#' towards infinity because the length of one degree of longitude tends towards
#' 0. For regions that span only a few degrees and are not too close to the
#' poles, setting the aspect ratio of the plot to the appropriate lat/lon ratio
#' approximates the usual mercator projection. This is what
#' `coord_quickmap()` does, and is much faster (particularly for complex
#' plots like [geom_tile()]) at the expense of correctness.
#'
#' @param projection projection to use, see
#' [mapproj::mapproject()] for list
#' @param ...,parameters Other arguments passed on to
#' [mapproj::mapproject()]. Use `...` for named parameters to
#' the projection, and `parameters` for unnamed parameters.
#' `...` is ignored if the `parameters` argument is present.
#' @param orientation projection orientation, which defaults to
#' `c(90, 0, mean(range(x)))`. This is not optimal for many
#' projections, so you will have to supply your own. See
#' [mapproj::mapproject()] for more information.
#' @param xlim,ylim Manually specific x/y limits (in degrees of
#' longitude/latitude)
#' @param clip Should drawing be clipped to the extent of the plot panel? A
#' setting of `"on"` (the default) means yes, and a setting of `"off"`
#' means no. For details, please see [`coord_cartesian()`].
#' @export
#' @examples
#' if (require("maps")) {
#' nz <- map_data("nz")
#' # Prepare a map of NZ
#' nzmap <- ggplot(nz, aes(x = long, y = lat, group = group)) +
#' geom_polygon(fill = "white", colour = "black")
#'
#' # Plot it in cartesian coordinates
#' nzmap
#' }
#'
#' if (require("maps")) {
#' # With correct mercator projection
#' nzmap + coord_map()
#' }
#'
#' if (require("maps")) {
#' # With the aspect ratio approximation
#' nzmap + coord_quickmap()
#' }
#'
#' if (require("maps")) {
#' # Other projections
#' nzmap + coord_map("azequalarea", orientation = c(-36.92, 174.6, 0))
#' }
#'
#' if (require("maps")) {
#' states <- map_data("state")
#' usamap <- ggplot(states, aes(long, lat, group = group)) +
#' geom_polygon(fill = "white", colour = "black")
#'
#' # Use cartesian coordinates
#' usamap
#' }
#'
#' if (require("maps")) {
#' # With mercator projection
#' usamap + coord_map()
#' }
#'
#' if (require("maps")) {
#' # See ?mapproject for coordinate systems and their parameters
#' usamap + coord_map("gilbert")
#' }
#'
#' if (require("maps")) {
#' # For most projections, you'll need to set the orientation yourself
#' # as the automatic selection done by mapproject is not available to
#' # ggplot
#' usamap + coord_map("orthographic")
#' }
#'
#' if (require("maps")) {
#' usamap + coord_map("conic", lat0 = 30)
#' }
#'
#' if (require("maps")) {
#' usamap + coord_map("bonne", lat0 = 50)
#' }
#'
#' \dontrun{
#' if (require("maps")) {
#' # World map, using geom_path instead of geom_polygon
#' world <- map_data("world")
#' worldmap <- ggplot(world, aes(x = long, y = lat, group = group)) +
#' geom_path() +
#' scale_y_continuous(breaks = (-2:2) * 30) +
#' scale_x_continuous(breaks = (-4:4) * 45)
#'
#' # Orthographic projection with default orientation (looking down at North pole)
#' worldmap + coord_map("ortho")
#' }
#'
#' if (require("maps")) {
#' # Looking up up at South Pole
#' worldmap + coord_map("ortho", orientation = c(-90, 0, 0))
#' }
#'
#' if (require("maps")) {
#' # Centered on New York (currently has issues with closing polygons)
#' worldmap + coord_map("ortho", orientation = c(41, -74, 0))
#' }
#' }
coord_map <- function(projection="mercator", ..., parameters = NULL, orientation = NULL, xlim = NULL, ylim = NULL, clip = "on") {
if (is.null(parameters)) {
params <- list(...)
} else {
params <- parameters
}
ggproto(NULL, CoordMap,
projection = projection,
orientation = orientation,
limits = list(x = xlim, y = ylim),
params = params,
clip = clip
)
}
#' @rdname ggplot2-ggproto
#' @format NULL
#' @usage NULL
#' @export
CoordMap <- ggproto("CoordMap", Coord,
transform = function(self, data, panel_params) {
trans <- mproject(self, data$x, data$y, panel_params$orientation)
out <- cunion(trans[c("x", "y")], data)
out$x <- rescale(out$x, 0:1, panel_params$x.proj)
out$y <- rescale(out$y, 0:1, panel_params$y.proj)
# mproject() converts Inf to NA, so we need to restore them from data.
out$x[is.infinite(data$x)] <- squish_infinite(data$x)
out$y[is.infinite(data$y)] <- squish_infinite(data$y)
out
},
backtransform_range = function(panel_params) {
# range is stored in data coordinates and doesn't have to be back-transformed
list(x = panel_params$x.range, y = panel_params$y.range)
},
range = function(panel_params) {
# Range in projected coordinates:
# list(x = panel_params$x.proj, y = panel_params$y.proj)
# However, coord_map() does never really work with transformed coordinates,
# so return unprojected data coordinates here
list(x = panel_params$x.range, y = panel_params$y.range)
},
distance = function(x, y, panel_params) {
max_dist <- dist_central_angle(panel_params$x.range, panel_params$y.range)
dist_central_angle(x, y) / max_dist
},
aspect = function(ranges) {
diff(ranges$y.proj) / diff(ranges$x.proj)
},
setup_panel_params = function(self, scale_x, scale_y, params = list()) {
# range in scale
ranges <- list()
for (n in c("x", "y")) {
scale <- get(paste0("scale_", n))
limits <- self$limits[[n]]
range <- expand_limits_scale(scale, default_expansion(scale), coord_limits = limits)
ranges[[n]] <- range
}
orientation <- self$orientation %||% c(90, 0, mean(ranges$x))
# Increase chances of creating valid boundary region
grid <- expand.grid(
x = seq(ranges$x[1], ranges$x[2], length.out = 50),
y = seq(ranges$y[1], ranges$y[2], length.out = 50)
)
ret <- list(x = list(), y = list())
# range in map
proj <- mproject(self, grid$x, grid$y, orientation)$range
ret$x$proj <- proj[1:2]
ret$y$proj <- proj[3:4]
for (n in c("x", "y")) {
out <- get(paste0("scale_", n))$break_info(ranges[[n]])
ret[[n]]$range <- out$range
ret[[n]]$major <- out$major_source
ret[[n]]$minor <- out$minor_source
ret[[n]]$labels <- out$labels
}
details <- list(
orientation = orientation,
x.range = ret$x$range, y.range = ret$y$range,
x.proj = ret$x$proj, y.proj = ret$y$proj,
x.major = ret$x$major, x.minor = ret$x$minor, x.labels = ret$x$labels,
y.major = ret$y$major, y.minor = ret$y$minor, y.labels = ret$y$labels,
x.arrange = scale_x$axis_order(), y.arrange = scale_y$axis_order()
)
details
},
render_bg = function(self, panel_params, theme) {
xrange <- expand_range(panel_params$x.range, 0.2)
yrange <- expand_range(panel_params$y.range, 0.2)
# Limit ranges so that lines don't wrap around globe
xmid <- mean(xrange)
ymid <- mean(yrange)
xrange[xrange < xmid - 180] <- xmid - 180
xrange[xrange > xmid + 180] <- xmid + 180
yrange[yrange < ymid - 90] <- ymid - 90
yrange[yrange > ymid + 90] <- ymid + 90
xgrid <- with(panel_params, expand.grid(
y = c(seq(yrange[1], yrange[2], length.out = 50), NA),
x = x.major
))
ygrid <- with(panel_params, expand.grid(
x = c(seq(xrange[1], xrange[2], length.out = 50), NA),
y = y.major
))
xlines <- self$transform(xgrid, panel_params)
ylines <- self$transform(ygrid, panel_params)
if (nrow(xlines) > 0) {
grob.xlines <- element_render(
theme, "panel.grid.major.x",
xlines$x, xlines$y, default.units = "native"
)
} else {
grob.xlines <- zeroGrob()
}
if (nrow(ylines) > 0) {
grob.ylines <- element_render(
theme, "panel.grid.major.y",
ylines$x, ylines$y, default.units = "native"
)
} else {
grob.ylines <- zeroGrob()
}
ggname("grill", grobTree(
element_render(theme, "panel.background"),
grob.xlines, grob.ylines
))
},
render_axis_h = function(self, panel_params, theme) {
arrange <- panel_params$x.arrange %||% c("secondary", "primary")
if (is.null(panel_params$x.major)) {
return(list(
top = zeroGrob(),
bottom = zeroGrob()
))
}
x_intercept <- with(panel_params, data_frame0(
x = x.major,
y = y.range[1],
.size = length(x.major)
))
pos <- self$transform(x_intercept, panel_params)
axes <- list(
top = draw_axis(pos$x, panel_params$x.labels, "top", theme),
bottom = draw_axis(pos$x, panel_params$x.labels, "bottom", theme)
)
axes[[which(arrange == "secondary")]] <- zeroGrob()
axes
},
render_axis_v = function(self, panel_params, theme) {
arrange <- panel_params$y.arrange %||% c("primary", "secondary")
if (is.null(panel_params$y.major)) {
return(list(
left = zeroGrob(),
right = zeroGrob()
))
}
x_intercept <- with(panel_params, data_frame0(
x = x.range[1],
y = y.major,
.size = length(y.major)
))
pos <- self$transform(x_intercept, panel_params)
axes <- list(
left = draw_axis(pos$y, panel_params$y.labels, "left", theme),
right = draw_axis(pos$y, panel_params$y.labels, "right", theme)
)
axes[[which(arrange == "secondary")]] <- zeroGrob()
axes
}
)
mproject <- function(coord, x, y, orientation) {
check_installed("mapproj", reason = "for `coord_map()`")
suppressWarnings(mapproj::mapproject(x, y,
projection = coord$projection,
parameters = coord$params,
orientation = orientation
))
}
|