Nothing
## mapping.R
##
## Copyright (C) 2015, Finn Lindgren
# old_globeproj ####
#' @title Old globe projection methods
#' @description Deprecated globe projection methods that may be removed in the
#' future
#' @keywords internal
#' @param type Projection type, see [.globeproj.types()]
#' @param orient long,lat,rotation
#' @param xlim x-axis limits
#' @param ylim y-axis limits
#' @param scale x- and y- scaling factors
#' @author Finn Lindgren
#' @name globeproj
#' @export
#' @examples
#' old_globeproj("mollweide")
#'
old_globeproj <- function(type = NULL,
orient = NULL,
xlim = NULL,
ylim = NULL,
scale = NULL) {
if (missing(type)) {
return(.globeproj.types())
}
type <- match.arg(type, .globeproj.types())
name <- type
if (identical(type, "lambert")) {
type <- "orthocyl"
if (is.null(scale)) scale <- c(1, 1)
if (is.null(orient)) orient <- c(0, 0, 0, 0)
} else if (identical(type, "gall-peters")) {
type <- "orthocyl"
if (is.null(scale)) scale <- c(1, 1)
scale <- scale * c(1 * cos(pi / 4), 1 / cos(pi / 4))
if (is.null(orient)) orient <- c(0, 0, 0, 0)
}
if (identical(type, "longlat")) {
if (is.null(orient)) orient <- c(0, 0, 0, 0)
if (is.null(scale)) scale <- c(1, 1)
if (is.null(xlim)) xlim <- c(-180, 180) * scale[1]
if (is.null(ylim)) ylim <- c(-90, 90) * scale[2]
} else if (identical(type, "orthocyl")) {
if (is.null(orient)) orient <- c(0, 0, 0, 0)
if (is.null(scale)) scale <- c(1, 1)
if (is.null(xlim)) xlim <- c(-pi, pi) * scale[1]
if (is.null(ylim)) ylim <- c(-1, 1) * scale[2]
} else if (identical(type, "mollweide")) {
if (is.null(orient)) orient <- c(0, 0, 0, 0)
if (is.null(scale)) scale <- c(1, 1) / sqrt(2)
if (is.null(xlim)) xlim <- c(-2, 2) * sqrt(2) * scale[1]
if (is.null(ylim)) ylim <- c(-1, 1) * sqrt(2) * scale[2]
} else if (identical(type, "hammer")) {
if (is.null(orient)) orient <- c(0, 0, 0, 0)
if (is.null(scale)) scale <- c(1, 1) / sqrt(2)
if (is.null(xlim)) xlim <- c(-2, 2) * sqrt(2) * scale[1]
if (is.null(ylim)) ylim <- c(-1, 1) * sqrt(2) * scale[2]
}
while (length(orient) < 4) {
orient <- c(orient, 0)
}
x <- structure(
list(
name = name, type = type,
orient = orient,
xlim = xlim, ylim = ylim, scale = scale
),
class = "globeproj"
)
.validobject_globeproj(x)
x
}
##' @describeIn globeproj Types of globe projections
##' @param x A [globeproj] object.
##' @param \dots Not used.
##' @return A vector of names of available projection types; "longlat",
##' "mollweide", "hammer", "orthocyl", "lambert", "gall-peters"
##' @author Finn Lindgren
.globeproj.types <-
function(x, ...) {
c(
"longlat",
"mollweide", "hammer",
"orthocyl",
"lambert", "gall-peters"
)
}
valid_globeproj <- function(object, ...) {
msg <- c()
if (!(object$type %in% .globeproj.types())) {
msg <- c(
msg,
paste("'type' should be one of (",
paste("'", .globeproj.types(), "'",
collapse = ", ", sep = ""
),
"), not '", object$type, "'.",
sep = ""
)
)
}
if (length(object$xlim) != 2) {
msg <- c(
msg,
paste("'xlim' should have length 2, not ",
length(object$xlim), ".",
sep = ""
)
)
}
if (length(object$ylim) != 2) {
msg <- c(
msg,
paste("'ylim' should have length 2, not ",
length(object$ylim), ".",
sep = ""
)
)
}
if (length(msg) == 0) {
TRUE
} else {
msg
}
}
.validobject_globeproj <- function(x, test = FALSE, ...) {
## Mostly from validObject
errors <- valid_globeproj(x)
if (length(errors) > 1 || is.character(errors[1])) {
if (test) {
errors
} else {
msg <- gettextf("invalid class %s object", dQuote(class(x)))
if (length(errors) > 1L) {
stop(paste(paste0(msg, ":"), paste(seq_along(errors),
errors,
sep = ": "
),
collapse = "\n"
), domain = NA)
} else {
stop(msg, ": ", errors, domain = NA)
}
}
} else {
TRUE
}
}
.clip_globeproj <- function(x, coords) {
## Clip and generate a list of Line objects
thelines <- list()
## Rudimentary cutting:
toolong <-
which(c(
TRUE,
diff(coords[, 1])^2 / diff(x$xlim)^2
+ diff(coords[, 2])^2 / diff(x$ylim)^2
> 0.1^2,
TRUE
))
start <- toolong[-length(toolong)]
ending <- toolong[-1] - 1
for (i in seq_along(start)) {
coords1 <- coords[start[i]:ending[i], 1]
coords2 <- coords[start[i]:ending[i], 2]
thelines <-
c(
thelines,
list(sp::Line(cbind(coords1, coords2)))
)
}
thelines
}
# plot_Polyset ####
#' @title Plot a projected PolySet
#' @param x A PolySet (see `PBSmapping`)
#' @param projection a `globeproj` objcet
#' @param add logical; If TRUE, add to existing plot
#' @param \dots Additional parameters passed on to `sp::plot`
#' @details This function is a legacy method for plotting data from the
#' `PBSmapping` package, with
#'
#' ```
#' data(worldLL, package = "PBSmapping")
#' plot_PolySet(worldLL, old_globeproj("longlat"), add = FALSE)
#' ```
#'
#' To avoid a dependency on the `PBSmapping` package, the example below
#' constructs a synthetic object of the same format.
#' @export
#' @seealso [old_globeproj()]
#' @examples
#' world_example <- data.frame(
#' PID = c(0L, 0L, 0L, 1L, 1L),
#' POS = c(1L, 2L, 3L, 1L, 2L),
#' X = c(10, 20, 30, 15, 25),
#' Y = c(10, 15, 70, -40, -50)
#' )
#' plot_PolySet(world_example, old_globeproj("longlat"), add = FALSE)
#' @keywords internal
#' @returns An (invisible) `sp` object of projected lines
plot_PolySet <- function(x, projection, add = FALSE, ...) {
coords <- old_project(projection, cbind(x$X, x$Y))
proj <-
sp::SpatialLines(list(sp::Lines(
unlist(
lapply(
unique(x$PID),
function(k) {
.clip_globeproj(
projection,
coords[sum(x$PID < k) +
x$POS[x$PID == k], , drop = FALSE]
)
}
),
recursive = FALSE
),
ID = "shape"
)))
if (add) {
sp::plot(proj, add = TRUE, ...)
} else {
sp::plot(proj, ...)
}
invisible(proj)
}
# old_outline ####
#' @param x A [globeproj] object
#' @param add logical; If TRUE, add to existing plot
#' @param do.plot logical; If try, do plotting
#' @param \dots Additional parameters passed to other methods
#' @return A
#' @author Finn Lindgren
#' @export old_outline
#' @aliases outline
#' @aliases old_outline
#' @rdname globeproj
old_outline <-
function(x, add = FALSE, do.plot = TRUE,
...) {
thebox <-
sp::SpatialPolygons(list(sp::Polygons(
list(sp::Polygon(
cbind(
c(
x$xlim[1], x$xlim[1],
x$xlim[2], x$xlim[2], x$xlim[1]
),
c(
x$ylim[1], x$ylim[2],
x$ylim[2], x$ylim[1], x$ylim[1]
)
)
)),
ID = "box"
)))
thebox <- sf::st_as_sfc(thebox)
if (x$type %in% c("longlat", "orthocyl")) {
theoutline <- thebox
} else if (x$type %in% c("mollweide", "hammer")) {
angle <- seq(0, 2 * pi, length = 361)
circle <- cbind(cos(angle), -sin(angle))
theboundary <-
sp::SpatialPolygons(list(sp::Polygons(
list(sp::Polygon(
cbind(
circle[, 1] * 2 * sqrt(2) * x$scale[1],
circle[, 2] * sqrt(2) * x$scale[2]
)
)),
ID = "boundary"
)))
theboundary <- sf::st_as_sfc(theboundary)
theoutline <- sf::st_intersection(theboundary, thebox)
}
theoutline <- sf::as_Spatial(theoutline)
if (do.plot) {
if (add) {
sp::plot(theoutline, add = TRUE, ...)
} else {
sp::plot(theoutline, ...)
}
}
invisible(theoutline)
}
# old_graticule ####
#' @rdname globeproj
#' @param x A [globeproj] object
#' @param n The number of graticules (n-long, n-lat) to compute
#' @param add logical; If TRUE, add to existing plot
#' @param do.plot logical; If TRUE, do plotting
#' @param \dots Additional parameters passed on to other methods
#' @return A
#' @author Finn Lindgren
#' @export old_graticule
#' @aliases graticule
#' @aliases old_graticule
old_graticule <-
function(x, n = c(24, 12), add = FALSE, do.plot = TRUE,
...) {
## Graticule
if (n[1] > 0) {
graticule1 <- floor(n[1] / 2)
lon <- ((-graticule1):graticule1) * 2 / n[1] * 180
lat <- seq(-90, 90, by = 2)
meridians <- expand.grid(lat, lon)[, 2:1]
## TODO: Add oblique support (requires cutting), etc
proj.mer.coords <- old_project(x, meridians)
proj.mer.coords1 <- matrix(
proj.mer.coords[, 1], length(lat),
length(lon)
)
proj.mer.coords2 <- matrix(
proj.mer.coords[, 2], length(lat),
length(lon)
)
proj.mer <-
sp::SpatialLines(list(sp::Lines(
unlist(
lapply(
seq_along(lon),
function(k) {
.clip_globeproj(x, cbind(
proj.mer.coords1[, k],
proj.mer.coords2[, k]
))
}
),
recursive = FALSE
),
ID = "meridians"
)))
if (do.plot) {
if (add) {
sp::plot(proj.mer, add = TRUE, ...)
} else {
sp::plot(proj.mer, ...)
add <- TRUE
}
}
} else {
proj.mer <- NULL
}
if (n[2] > 1) {
lon <- seq(-180, 180, by = 2)
lat <- seq(-90, 90, length = 1 + n[2])[-c(1, 1 + n[2])]
parallels <- expand.grid(lon, lat)
proj.par.coords <- old_project(x, parallels)
proj.par.coords1 <- matrix(
proj.par.coords[, 1], length(lon),
length(lat)
)
proj.par.coords2 <- matrix(
proj.par.coords[, 2], length(lon),
length(lat)
)
proj.par <-
sp::SpatialLines(list(sp::Lines(
unlist(
lapply(
seq_along(lat),
function(k) {
.clip_globeproj(x, cbind(
proj.par.coords1[, k],
proj.par.coords2[, k]
))
}
),
recursive = FALSE
),
ID = "parallels"
)))
if (do.plot) {
if (add) {
sp::plot(proj.par, add = TRUE, ...)
} else {
sp::plot(proj.par, ...)
}
}
} else {
proj.par <- NULL
}
invisible(list(meridians = proj.mer, parallels = proj.par))
}
# old_tissot ####
##' @rdname globeproj
##' @param x A [globeproj] object
##' @param n The number of Tissot indicatrices (n-long, n-lat) to compute
##' @param add logical; If TRUE, add to existing plot
##' @param do.plot logical; If TRUE, do plotting
##' @param \dots Additional parameters passed on to other methods
##' @return A
##' @author Finn Lindgren
##' @export old_tissot
##' @aliases tissot
##' @aliases old_tissot
old_tissot <-
function(x, n = c(12, 6), add = FALSE, do.plot = TRUE,
...) {
## Tissot
}
# plot_globeproj ####
#' @title Plot a globeproj object
#' @param x A [globeproj] object
#' @param xlim,ylim The x- and y-axis limits
#' @param outline logical
#' @param graticule The number of graticules (n-long, n-lat) to compute
#' @param tissot The number of Tissot indicatrices (n-long, n-lat) to compute
#' @param asp the aspect ratio. Default = 1
#' @param add logical; If `TRUE`, add to existing plot. Default: `FALSE`
#' @param \dots Additional parameters passed on to other methods
#' @return Nothing
#' @author Finn Lindgren
#' @export
#' @examples
#' proj <- old_globeproj("moll", orient = c(0, 0, 45))
#' plot_globeproj(proj, graticule = c(24, 12), add = FALSE, asp = 1, lty = 2, lwd = 0.5)
plot_globeproj <-
function(x, xlim = NULL, ylim = NULL,
outline = TRUE,
graticule = c(24, 12),
tissot = c(12, 6),
asp = 1,
add = FALSE,
...) {
if (is.null(xlim)) xlim <- x$xlim
if (is.null(ylim)) ylim <- x$ylim
if (!add) {
plot.new()
plot(NA, type = "n", xlim = xlim, ylim = ylim, asp = asp, ...)
}
## Outline
if (outline) {
old_outline(x, add = TRUE, ...)
}
## Graticule
old_graticule(x, n = graticule, add = TRUE, ...)
## Tissot
old_tissot(x, n = tissot, add = TRUE, asp = asp, ...)
invisible(NULL)
}
# old_limits ####
##' @describeIn globeproj Calculates projection axis limits
##'
##' @param x A [globeproj] object
##' @param loc Coordinates to be mapped.
##' @param \dots Additional parameters passed on to other methods
##' @return A list:
##' \item{xlim }{X axis limits in the map domain}
##' \item{ylim }{Y axis limits in the map domain}
##' @author Finn Lindgren
##' @export old_limits
##' @aliases limits
##' @aliases old_limits
old_limits <-
function(x,
loc = NULL,
...) {
if (!is(x, "globeproj")) {
stop("'x' must be a 'globeproj'")
}
lim <- list(xlim = x$xlim, ylim = x$ylim)
if (!is.null(loc)) {
locproj <- old_project(x, loc)
xlim <- range(locproj[, 1], na.rm = TRUE)
ylim <- range(locproj[, 2], na.rm = TRUE)
lim$xlim <- c(max(lim$xlim[1], xlim[1]), min(lim$xlim[2], xlim[2]))
lim$ylim <- c(max(lim$ylim[1], ylim[1]), min(lim$ylim[2], ylim[2]))
}
lim
}
# old_project ####
#' @rdname globeproj
#' @param x A [globeproj] object
#' @param loc Coordinates to be mapped.
#' @param inverse logical
#' @param \dots Additional parameters passed on to other methods
#' @return B
#' @author Finn Lindgren
#'
#' @export old_project
#' @aliases project
#' @aliases old_project
old_project <-
function(x, loc, inverse = FALSE, ...) {
if (!is(x, "globeproj")) {
stop("'x' must be a 'globeproj'")
}
if (inverse) {
## Convert coordinates to standardised scale
loc <- loc %*% diag(1 / x$scale, 2)
} else {
if (ncol(loc) == 2) {
## Standardise (long,lat) to Cartesian coordinates.
loc <-
cbind(
x = cos(loc[, 1] * pi / 180) * cos(loc[, 2] * pi / 180),
y = sin(loc[, 1] * pi / 180) * cos(loc[, 2] * pi / 180),
z = sin(loc[, 2] * pi / 180)
)
}
## Transform to oblique orientation
## 1) Rotate -orient[1] around (0,0,1)
## 2) Rotate +orient[2] around (0,1,0)
## 3) Rotate -orient[3] around (1,0,0)
## 3) Rotate -orient[4] around (0,0,1)
loc <- loc %*% fm_rotmat3213(c(-1, 1, -1, -1) * x$orient * pi / 180)
}
if (identical(x$type, "longlat")) {
if (inverse) {
proj <-
cbind(
x = cos(loc[, 1] * pi / 180) * cos(loc[, 2] * pi / 180),
y = sin(loc[, 1] * pi / 180) * cos(loc[, 2] * pi / 180),
z = sin(loc[, 2] * pi / 180)
)
} else {
proj <-
cbind(
x = atan2(loc[, 2], loc[, 1]) * 180 / pi,
y = asin(pmax(-1, pmin(+1, loc[, 3]))) * 180 / pi
)
}
} else if (identical(x$type, "orthocyl")) {
if (inverse) {
coslat <- sqrt(pmax(0, 1 - loc[, 2]^2))
proj <-
cbind(
x = cos(loc[, 1] * pi / 180) * coslat,
y = sin(loc[, 1] * pi / 180) * coslat,
z = loc[, 2]
)
} else {
proj <-
cbind(
x = atan2(loc[, 2], loc[, 1]),
y = loc[, 3]
)
}
} else if (identical(x$type, "mollweide")) {
if (inverse) {
ok <- ((loc[, 1]^2 / 2 + 2 * loc[, 2]^2) <= 4)
cos.theta <- sqrt(pmax(0, 1 - loc[ok, 2]^2 / 2))
theta <- atan2(loc[ok, 2] / sqrt(2), cos.theta)
sin.lat <- (2 * theta + sin(2 * theta)) / pi
cos.lat <- sqrt(pmax(0, 1 - sin.lat^2))
lon <- loc[ok, 1] / sqrt(2) * pi / 2 / (cos.theta + (cos.theta == 0))
lon[cos.theta == 0] <- pi / 2 * sign(theta[cos.theta == 0])
proj <- matrix(NA, nrow(loc), 3)
proj[ok, ] <- cbind(
x = cos(lon) * cos.lat,
y = sin(lon) * cos.lat,
z = sin.lat
)
} else {
lon <- atan2(loc[, 2], loc[, 1])
z <- pmin(1, pmax(-1, loc[, 3]))
sin.theta <- z
cos.theta <- sqrt(pmax(0, 1 - sin.theta^2))
## NR-solver for sin.theta.
## Typically finishes after at most 7 iterations.
## When cos.theta=0, sin.theta is already correct, +/- 1.
notok <- (cos.theta > 0)
for (k in 1:20) {
if (!any(notok)) {
break
}
delta <-
(atan2(sin.theta[notok], cos.theta[notok]) +
sin.theta[notok] * cos.theta[notok] - pi / 2 * z[notok]) /
(2 * cos.theta[notok])
sin.theta[notok] <- sin.theta[notok] - delta
cos.theta[notok] <- sqrt(1 - sin.theta[notok]^2)
notok[notok] <- (abs(delta) > 1e-14)
}
proj <- cbind(
x = 2 * lon / pi * cos.theta * sqrt(2),
y = sin.theta[, drop = TRUE] * sqrt(2)
)
}
} else if (identical(x$type, "hammer")) {
if (inverse) {
ok <- ((loc[, 1]^2 / 2 + 2 * loc[, 2]^2) <= 4)
z2 <- 1 - (loc[ok, 1] / 4)^2 - (loc[ok, 2] / 2)^2
z <- sqrt(z2)
lon <- 2 * atan(z * loc[ok, 1] / (2 * (2 * z2 - 1)))
sin.lat <- z * loc[ok, 2]
cos.lat <- sqrt(pmax(0, 1 - sin.lat^2))
proj <- matrix(NA, nrow(loc), 3)
proj[ok, ] <- cbind(
x = cos(lon) * cos.lat,
y = sin(lon) * cos.lat,
z = sin.lat
)
} else {
lon <- atan2(loc[, 2], loc[, 1])
sin.lat <- pmin(1, pmax(-1, loc[, 3]))
cos.lat <- sqrt(pmax(0, 1 - sin.lat^2))
cos.lon2 <- cos(lon / 2)
sin.lon2 <- sin(lon / 2)
scale <- sqrt(2) / sqrt(1 + cos.lat * cos.lon2)
proj <- cbind(
x = 2 * cos.lat * sin.lon2 * scale,
y = sin.lat * scale
)
}
}
if (inverse) {
## Transform back from oblique orientation
## 1) Rotate +orient[4] around (0,0,1)
## 2) Rotate +orient[3] around (1,0,0)
## 3) Rotate -orient[2] around (0,1,0)
## 4) Rotate +orient[1] around (0,0,1)
proj <- proj %*% fm_rotmat3123(c(1, -1, 1, 1) * x$orient * pi / 180)
} else {
proj <- cbind(x = proj[, 1] * x$scale[1], y = proj[, 2] * x$scale[2])
}
proj
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.