R/bb.R

In mtennekes/oldtmaptools: Old Thematic Map Tools

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

    ## 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
		} else current.projection <- get_proj4(current.projection, output = "crs")
		projection <- get_proj4(projection, output = "crs")

		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: try to cut the bounding box, such that it is feasible (i.e. corresponding to lon between -180 and 180 and lat between -90 and 90)
	    earth_end <- suppressWarnings(bb_earth(projection=current.projection))

	    if (is.null(earth_end)) {
	        sf_poly2 <- sf_poly
	    } else {
	        sf_poly2 <- tryCatch({
	            suppressMessages(sf::st_intersection(sf_poly, earth_end))
	        }, error=function(e){
	            sf_poly
	        })
	        if (is.null(sf_poly2) || !inherits(sf_poly2, c("sf", "sfc")) || length(sf_poly2) == 0) sf_poly2 <- sf_poly
	    }

	    # STEP 2: Extract the bounding box corner points and add intermediate points, which can be needed since the exterme values may not be corner points once they are projected. Create a SpatialPoints objects from these points.
	    co <- sf::st_coordinates(sf_poly2)[,1:2]
	    co2 <- apply(co, 2, function(v) {
	        n <- length(v)
	        c(v[1], rep(v[-n], each=4) + as.vector(sapply((v[-1] - v[-n]) / 4, function(w)cumsum(rep(w,4)))))
	    })

	    sf_pnts <- sf::st_sfc(sf::st_multipoint(co2), crs=current.projection)


	    # STEP 3: Reproject SpatialPoints object
	    tryCatch({
	        sf_pnts2_prj <- st_transform2(sf_pnts, crs=projection)
	    }, error=function(e) {
	        stop("Something went wrong with the bounding box. Please check the projection.", call.=FALSE)
	    })

	    # STEP 4: Get bounding box of reprojected object



	    b <- sf::st_bbox(sf_pnts2_prj)
	    is_prj <- is_projected(projection)
	} else {
	    is_prj <- if (is.na(current.projection)) {
	        !maybe_longlat(b)
	    } else is_projected(current.projection)

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



	## check if long lat coordinates are valid
	if (!is_prj) {
	    b[1:2] <- pmax(b[1:2], c(-180, -90))
	    b[3:4] <- pmin(b[3:4], c(180, 90))
	}

	output <- match.arg(output)

	if (output == "bbox") {
	    if (!inherits(b, "bbox")) {
	        b <- unname(b)
	        b <- st_bbox(c(xmin = b[1], ymin = b[2], xmax = b[3], ymax = b[4]), crs = st_crs(current.projection))
	    }
	} else if (output == "matrix") {
	    b <- matrix(b, ncol=2, dimnames = list(c("x", "y"), c("min", "max")))
	} else {
	    b <- extent(b[c(1,3,2,4)])
	}

	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))
        st_bbox(c(xmin = bb[1], ymin = bb[2], xmax = bb[3], ymax = bb[4]), crs = st_crs(NA))
    } else if (inherits(bb, "Extent")) {
        bb <- unname(as.vector(bb))
        st_bbox(c(xmin = bb[1], ymin = bb[3], xmax = bb[2], ymax = bb[4]), crs = 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 <- get_projection(x, output = "crs")
    } else if (inherits(x, "Spatial")) {
        b <- sfbb(attr(x, "bbox"))
        current.projection <- get_projection(x, output = "crs")
    } else if (inherits(x, c("sf", "sfc"))) {
        b <- st_bbox(x)
        current.projection <- 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 <- st_bbox(c(xmin = x[1], ymin = x[2], xmax = x[3], ymax = x[4]), crs = 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 <- st_bbox(c(xmin = xlim[1], ymin = ylim[1], xmax = xlim[2], ymax = ylim[2]), crs = st_crs(NA))

    }
    if (is.null(current.projection)) current.projection <- st_crs(NA)
    if (!is.na(current.projection)) {
        attr(b, "crs") <- current.projection
    } else if (!is.na(attr(b, "crs"))) {
        current.projection <- attr(b, "crs")
    }

    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
}