#' Overlay Multi-Polygon Objects
#'
#' Calculate the intersection or difference between two multi-polygon objects.
#'
#' @param x 'SpatialPolygons*'.
#' Multi-polygon object
#' @param y 'SpatialPolygons*' or 'Extent'.
#' Multi-polygon object
#' @param cmd 'character' string.
#' Specifying "\link{gIntersection}", the default, cuts out portions of the
#' \code{x} polygons that overlay the \code{y} polygons.
#' If "\link{gDifference}" is specified, only those portions of the \code{x} polygons
#' falling outside the \code{y} polygons are copied to the output polygons.
#' @param buffer.width 'numeric' number.
#' Expands or contracts the geometry of \code{y} to include the area within
#' the specified width, see \code{gBuffer}.
#' Specifying \code{NA}, the default, indicates no buffer.
#'
#' @details This function tests if the resulting geometry is valid.
#' If invalid, an attempt is made to make the geometry valid by zero-width buffering.
#'
#' @return An object of class 'SpatialPolygons*'.
#'
#' @author J.C. Fisher, U.S. Geological Survey, Idaho Water Science Center
#'
#' @keywords utilities
#'
#' @export
#'
#' @examples
#' m1a <- rbind(c(17.5, 55.1),
#' c(24.7, 55.0),
#' c(22.6, 61.1),
#' c(16.5, 59.7),
#' c(17.5, 55.1))
#' m1b <- m1a
#' m1b[, 1] <- m1b[, 1] + 11
#' p1 <- list(sp::Polygon(m1a, FALSE), sp::Polygon(m1b, FALSE))
#' p1 <- sp::SpatialPolygons(list(sp::Polygons(p1, 1)))
#' sp::plot(p1, col = "blue")
#'
#' m2a <- rbind(c(19.6, 60.0),
#' c(35.7, 58.8),
#' c(28.2, 64.4),
#' c(19.6, 60.0))
#' m2b <- rbind(c(20.6, 56.2),
#' c(30.9, 53.8),
#' c(27.3, 51.4),
#' c(20.6, 56.2))
#' p2 <- list(sp::Polygon(m2a, FALSE), sp::Polygon(m2b, FALSE))
#' p2 <- sp::SpatialPolygons(list(sp::Polygons(p2, 2)))
#' sp::plot(p2, col = "red", add = TRUE)
#'
#' p <- SetPolygons(p1, p2, "gIntersection")
#' sp::plot(p, col = "green", add = TRUE)
#'
#' p <- SetPolygons(p2, p1, "gDifference")
#' sp::plot(p, col = "purple", add = TRUE)
#'
SetPolygons <- function(x, y, cmd=c("gIntersection", "gDifference"), buffer.width=NA) {
stopifnot(inherits(x, c("SpatialPolygons", "SpatialPolygonsDataFrame")))
stopifnot(inherits(y, c("SpatialPolygons", "SpatialPolygonsDataFrame", "Extent")))
cmd <- match.arg(cmd)
checkmate::assertNumber(buffer.width, na.ok=TRUE, finite=TRUE)
if (inherits(y, "Extent")) {
crds <- cbind(c(y[1:2], y[2:1], y[1]), c(rep(y[3], 2), rep(y[4], 2), y[3]))
y <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(crds)), "bbox")),
proj4string=raster::crs(x))
}
if (inherits(x, "SpatialPolygonsDataFrame")) {
d <- x@data
rownames(d) <- sapply(methods::slot(x, "polygons"), function(i) methods::slot(i, "ID"))
} else {
d <- NULL
}
x <- methods::as(x, "SpatialPolygons")
y <- methods::as(y, "SpatialPolygons")
y <- y[which(apply(rgeos::gIntersects(y, x, byid=TRUE), 2, any)), ]
is_intersecting <- rgeos::gIntersects(x, y, byid=TRUE)
z <- lapply(seq_along(x), function (i) {
if (any(is_intersecting[, i])) {
y_intersect <- y[is_intersecting[, i]]
if (is.numeric(buffer.width))
y_intersect <- rgeos::gBuffer(y_intersect, width=buffer.width)
suppressMessages(suppressWarnings({
spgeom2 <- rgeos::gUnaryUnion(y_intersect, checkValidity=2L)
if (cmd == "gIntersection")
x_geo <- rgeos::gIntersection(x[i], spgeom2, byid=TRUE, checkValidity=2L)
else
x_geo <- rgeos::gDifference(x[i], spgeom2, byid=TRUE, checkValidity=2L)
if (inherits(x_geo, "SpatialCollections"))
x_geo <- rgeos::gUnaryUnion(x_geo@polyobj, checkValidity=2L)
}))
if (is.null(x_geo)) return(NULL)
p <- x_geo@polygons[[1]]
methods::slot(p, "ID") <- methods::slot(x[i]@polygons[[1]], "ID")
} else {
p <- if (cmd == "gIntersection") NULL else x[i]@polygons[[1]]
}
p
})
is_retained <- !vapply(z, is.null, TRUE)
z <- sp::SpatialPolygons(z[is_retained], proj4string=raster::crs(x))
if (inherits(d, "data.frame")) {
d <- d[is_retained, , drop=FALSE]
z <- sp::SpatialPolygonsDataFrame(z, d, match.ID=TRUE)
}
z
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.