R/bb.R

In tmaptools: Thematic Map Tools

#' 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}} and \url{https://wiki.openstreetmap.org/wiki/Nominatim}.}
#' 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.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.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.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)
}