#' Bounding box generator

#'

#' Swiss army knife for bounding boxes. Modify an existing bounding box or create a new bounding box from scratch. See details.

#'

#' An existing bounding box (defined by \code{x}) can be modified as follows:

#' \itemize{

#' \item Using the extension factor \code{ext}.

#' \item Changing the width and height with \code{width} and \code{height}. The argument \code{relavitve} determines whether relative or absolute values are used.

#' \item Setting the x and y limits. The argument \code{relavitve} determines whether relative or absolute values are used.}

#'

#' A new bounding box can be created from scratch as follows:

#' \itemize{

#' \item Using the extension factor \code{ext}.

#' \item Setting the center coorinates \code{cx} and \code{cy}, together with the \code{width} and \code{height}.

#' \item Setting the x and y limits \code{xlim} and \code{ylim}

#' }

#'

#' @param x One of the following:

#' \itemize{

#' \item A shape from class \code{\link[sf:sf]{sf}}, \code{\link[stars:st_as_stars]{stars}}, \code{sp}, or \code{raster}.

#' \item A bounding box (\code{\link[sf:st_bbox]{st_bbox}}, \code{Extent} (\code{raster} package, which will no longer be supported in the future versions), numeric vector of 4 (default order: xmin, ymin, xmax, ymax), or a 2x2 matrix).

#' \item Open Street Map search query. The bounding is automatically generated by querying \code{x} from Open Street Map Nominatim. See \code{\link{geocode_OSM}}}.

# If \code{x} is not specified, a bounding box can be created from scratch (see details).

#' @param ext Extension factor of the bounding box. If 1, the bounding box is unchanged. Values smaller than 1 reduces the bounding box, and values larger than 1 enlarges the bounding box. This argument is a shortcut for both \code{width} and \code{height} with \code{relative=TRUE}. If a negative value is specified, then the shortest side of the bounding box (so width or height) is extended with \code{ext}, and the longest side is extended with the same absolute value. This is especially useful for bounding boxes with very low or high aspect ratios.

#' @param cx center x coordinate

#' @param cy center y coordinate

#' @param width width of the bounding box. These are either absolute or relative (depending on the argument \code{relative}).

#' @param height height of the bounding box. These are either absolute or relative (depending on the argument \code{relative}).

#' @param xlim limits of the x-axis. These are either absolute or relative (depending on the argument \code{relative}).

#' @param ylim limits of the y-axis. See \code{xlim}.

#' @param relative boolean that determines whether relative values are used for \code{width}, \code{height}, \code{xlim} and \code{ylim} or absolute. If \code{x} is unspecified, \code{relative} is set to \code{"FALSE"}.

#' @param asp.target target aspect ratio, which is width/height, of the returned bounding box.

#' @param asp.limit maximum aspect ratio, which is width/height. Number greater than or equal to 1. For landscape bounding boxes, \code{1/asp.limit} will be used. The returned bounding box will have an aspect ratio between \code{1/asp.limit} and \code{asp.limit}.

#' @param current.projection projection that corresponds to the bounding box specified by \code{x}.

#' @param projection projection to transform the bounding box to.

#' @param output output format of the bounding box, one of:

#' \itemize{

#' \item \code{"bbox"} a \code{sf::bbox} object, which is a numeric vector of 4: xmin, ymin, xmax, ymax. This representation used by the \code{\link[sf:sf]{sf}} package.

#' \item \code{"matrix"} a 2 by 2 numeric matrix, where the rows correspond to x and y, and the columns to min and max. This representation used by the \code{sp} package.

#' \item \code{"extent"} an \code{raster::extent} object, which is a numeric vector of 4: xmin, xmax, ymin, ymax. This representation used by the \code{raster} package.

#' }

#' @return bounding box (see argument \code{output})

#' @import sf

#' @importFrom XML xmlTreeParse xmlChildren xmlRoot xmlAttrs

#' @example ./examples/bb.R

#' @seealso \code{\link{geocode_OSM}}

#' @export

bb <- function(x=NA, ext=NULL, cx=NULL, cy=NULL, width=NULL, height=NULL, xlim=NULL, ylim=NULL, relative = FALSE, asp.target = NULL, asp.limit = NULL, current.projection=NULL, projection=NULL, output = c("bbox", "matrix", "extent")) {

# check projections

if (!is.null(current.projection)) current.projection <- sf::st_crs(current.projection)

if (!is.null(projection)) projection <- sf::st_crs(projection)

## get unprocessed bounding box

res <- get_bb(x, cx=cx, cy=cy, width=width, height=height, xlim=xlim, ylim=ylim, current.projection=current.projection)

b <- res$b

cx <- res$cx

cy <- res$cy

width <- res$width

height <- res$height

xlim <- res$xlim

ylim <- res$ylim

current.projection <- res$current.projection

## impute cx and cy

if (!is.character(x)) {

if (is.null(cx)) cx <- mean(b[c(1,3)])

if (is.null(cy)) cy <- mean(b[c(2,4)])

}

## translate ext to width and height

steps <- b[3:4] - b[1:2]

if (!missing(ext)) {

relative <- TRUE

if (ext > 0) {

width <- ext

height <- ext

} else {

if (steps[1] > steps[2]) {

height <- -ext

fact <- (steps[2]/ steps[1])

if (is.nan(fact)) fact <- 1

width <- 1 + (-ext-1) * fact

} else {

width <- -ext

fact <- (steps[1]/ steps[2])

if (is.nan(fact)) fact <- 1

height <- 1 + (-ext-1) * fact

}

}

}

## modify bb

if (relative) {

xlim <- if (!is.null(xlim)) {

b[1] + xlim * steps[1]

} else if (!is.null(width)) {

c(cx - (width/2) * steps[1],

cx + (width/2) * steps[1])

} else {

b[c(1,3)]

}

ylim <- if (!is.null(ylim)) {

b[2] + ylim * steps[2]

} else if (!is.null(height)) {

c(cy - (height/2) * steps[2],

cy + (height/2) * steps[2])

} else {

b[c(2,4)]

}

} else {

if (!is.null(width)) {

xlim <- c(cx - (width/2),

cx + (width/2))

} else if (is.null(xlim)) {

xlim <- b[c(1,3)]

}

if (!is.null(height)) {

ylim <- c(cy - (height/2),

cy + (height/2))

} else if (is.null(ylim)) {

ylim <- b[c(2,4)]

}

}

## create bb

b <- c(xlim[1], ylim[1], xlim[2], ylim[2])

## reproject bb

if (!missing(projection)) {

#####################################################

# Reproject bounding box

#####################################################

if (is.na(current.projection)) {

if (!maybe_longlat(b)) {

stop("Current projection unknown. Please specify the projection.")

}

warning("Current projection unknown. Long lat coordinates (wgs84) assumed.", call. = FALSE)

current.projection <- .crs_longlat

}

mx = mean(b[c(1,3)])

my = mean(b[c(2,4)])

ls1 = sf::st_linestring(rbind(c(b[1], b[2]), c(b[3], b[2]),

c(b[3], b[4]), c(b[1], b[4]), c(b[1], b[2])))

ls2 = sf::st_linestring(rbind(c(b[1], b[2]), c(b[3], b[4]),

c(b[1], b[4]), c(b[3], b[2]), c(b[1], b[2])))

ls3 = sf::st_linestring(rbind(c(mx, b[2]), c(mx, b[4]),

c(b[1], my), c(b[3], my), c(mx, b[4]), c(b[3], my), c(mx, b[2]), c(b[1], my)))

sf_lns = sf::st_sfc(ls1, ls2,ls3)

sf_lns = sf::st_segmentize(sf_lns, sf::st_length(sf_lns)[1]/100)

sf::st_crs(sf_lns) = current.projection

#sf_poly <- sf::st_sfc(sf::st_polygon(list(matrix(c(b[1], b[2], b[1], b[4], b[3], b[4], b[3], b[2], b[1], b[2]), byrow = TRUE, ncol = 2))), crs=current.projection)

# STEP 1: check if poly can be reprojected, and if not, cut the bounding box such that lon between -180 and 180 and lat between -90 and 90

#sf_poly2 <- sf_poly

sf_lns_prj <- sf::st_transform(sf_lns, crs = projection, partial = TRUE)

is_prj <- !sf::st_is_longlat(projection)

#b <- sf::st_bbox(sf_lns_prj[!sf::st_is_empty(sf_lns_prj),])

b <- sf::st_bbox(sf_lns_prj)

if (any(!is.finite(b))) {

if (is_prj) stop("Unable to determine bounding for projected crs.", call. = FALSE)

b[1:4] <- c(-180, -90, 180, 90)

}

} else {

is_prj <- if (is.na(current.projection)) {

!maybe_longlat(b)

} else !sf::st_is_longlat(current.projection)

if (is.na(current.projection) && !is_prj) current.projection <- .crs_longlat

}

## check if lat coordinates are valid (long values may exceed the datetime boundary)

if (!is_prj) {

b[2] <- pmax(b[2], -90)

b[4] <- pmin(b[4], 90)

# b[1:2] <- pmax(b[1:2], c(-180, -90))

# b[3:4] <- pmin(b[3:4], c(180, 90))

}

## limit aspect ratio

if (!is.null(asp.target)) {

bbh <- unname(b[4] - b[2])

bbw <- unname(b[3] - b[1])

if (asp.target < (bbw / bbh)) {

# increase height

cy <- (b[2] + b[4]) / 2

b[2] <- cy - bbw / (asp.target * 2)

b[4] <- cy + bbw / (asp.target * 2)

} else if (asp.target > (bbw / bbh)) {

# increase width

cx <- (b[1] + b[3]) / 2

b[1] <- cx - bbh * (asp.target / 2)

b[3] <- cx + bbh * (asp.target / 2)

}

} else if (!is.null(asp.limit)) {

bbh <- unname(b[4] - b[2])

bbw <- unname(b[3] - b[1])

if (bbh < 1e-8 || bbh < (bbw / asp.limit)) {

# increase height such that aspect ratio is at most 10. Use lowerbound of 1 for height.

cy <- (b[2] + b[4]) / 2

b[2] <- cy - max(bbw / (asp.limit * 2), .5)

b[4] <- cy + max(bbw / (asp.limit * 2), .5)

}

if (bbw < 1e-8 || bbw < (bbh / asp.limit)) {

# increase width such that aspect ratio is at least 1/10. Use lowerbound of 1 for width.

cx <- (b[1] + b[3]) / 2

b[1] <- cx - max(bbh / (asp.limit * 2), .5)

b[3] <- cx + max(bbh / (asp.limit * 2), .5)

}

}

output <- match.arg(output)

if (output == "bbox") {

if (!inherits(b, "bbox")) {

b <- unname(b)

b <- sf::st_bbox(c(xmin = b[1], ymin = b[2], xmax = b[3], ymax = b[4]), crs = sf::st_crs(current.projection))

}

} else if (output == "matrix") {

b <- matrix(b, ncol=2, dimnames = list(c("x", "y"), c("min", "max")))

} else {

if (requireNamespace("raster")) {

warning("raster package functionality in tmaptools will no longer be supported from the next version onwards")

b <- raster::extent(b[c(1,3,2,4)])

} else stop("raster package needed")

}

b

}

get_sf_bbox <- function(shp) {

matrix(attr(shp[[attr(shp, "sf_column")]], "bbox"), ncol=2, dimnames = list(c("x", "y"), c("min", "max")))

}

sfbb <- function(bb) {

if (is.matrix(bb)) {

bb <- unname(as.vector(bb))

sf::st_bbox(c(xmin = bb[1], ymin = bb[2], xmax = bb[3], ymax = bb[4]), crs = sf::st_crs(NA))

} else if (inherits(bb, "Extent")) {

bb <- unname(as.vector(bb))

sf::st_bbox(c(xmin = bb[1], ymin = bb[3], xmax = bb[2], ymax = bb[4]), crs = sf::st_crs(NA))

} else stop("bb not 2x2 matrix nor extent object")

}

get_bb <- function(x, cx=NULL, cy=NULL, width=NULL, height=NULL, xlim=NULL, ylim=NULL, current.projection=NULL) {

if (is.character(x)) {

res <- as.vector(geocode_OSM(x))

b <- res$bbox

cx <- res$coords[1]

cy <- res$coords[2]

current.projection <- .crs_longlat

} else if (inherits(x, "Extent")) {

b <- sfbb(x)

} else if (inherits(x, "Raster")) {

b <- sfbb(attr(x, "extent"))

current.projection <- sf::st_crs(x)

} else if (inherits(x, "Spatial")) {

b <- sfbb(attr(x, "bbox"))

current.projection <- sf::st_crs(x)

} else if (inherits(x, c("sf", "sfc", "stars"))) {

b <- sf::st_bbox(x)

current.projection <- sf::st_crs(x)

} else if (is.matrix(x) && length(x)==4) {

b <- sfbb(x)

} else if (inherits(x, "bbox")) {

b <- x

} else if (is.vector(x) && length(x)==4) {

x <- unname(check_bb_order(x))

b <- sf::st_bbox(c(xmin = x[1], ymin = x[2], xmax = x[3], ymax = x[4]), crs = sf::st_crs(NA))

} else if (!is.na(x)[1]) {

stop("Incorrect x argument")

} else {

if ((is.null(xlim) && (is.null(width) || is.null(cx))) || (is.null(xlim) && (is.null(height) || is.null(cy))))

stop("Argument x is missing. Please specify x, or {xlim and ylim}, or {width, height, cx, and cy}.")

## create new bounding box

if (is.null(xlim)) xlim <- cx + c(-.5, .5) * width

if (is.null(ylim)) ylim <- cy + c(-.5, .5) * height

b <- sf::st_bbox(c(xmin = xlim[1], ymin = ylim[1], xmax = xlim[2], ymax = ylim[2]), crs = sf::st_crs(NA))

}

if (is.null(current.projection)) current.projection <- sf::st_crs(NA)

if (!is.na(current.projection)) {

attr(b, "crs") <- current.projection

} else if (!is.na(sf::st_crs(b))) {

current.projection <- sf::st_crs(b)

}

list(b=b,

cx=cx,

cy=cy,

width=width,

height=height,

xlim=xlim,

ylim=ylim,

current.projection=current.projection)

}

check_bb_order <- function(x) {

if (((x[1] > x[3]) || (x[2] > x[4])) && x[1] < x[2] && x[3] < x[4]) {

message("Bounding box format automatically changed from [xmin, xmax, ymin, ymax] to [xmin, ymin, xmax, ymax]")

x[c(1,3,2,4)]

} else x

}

maybe_longlat <- function(bb) {

(bb[1] >= -180.1 && bb[3] <= 180.1 && bb[2] >= -90.1 && bb[4] <= 90.1)

}