Nothing
#' Convert Spatial Grids to Polygons
#'
#' Convert gridded spatial data to spatial polygons.
#' Image files created with spatial polygons are reduced in size,
#' can easily be transformed from one coordinate reference system to another,
#' and result in much "cleaner" images when plotted.
#'
#' @param grd 'SpatialGridDataFrame', 'SpatialPixelsDataFrame', or 'Raster*'.
#' Spatial grid
#' @param zcol 'character' string or 'integer' count.
#' Layer to extract from a multi-layer spatial grid.
#' @param level 'logical' flag.
#' If true, a set of levels is used to partition the range of attribute values,
#' its default is false.
#' @param at 'numeric' vector.
#' Breakpoints along the range of attribute values.
#' @param cuts 'integer' count.
#' Number of levels the range of attribute values would be divided into.
#' @param pretty 'logical' flag.
#' Whether to use pretty cut locations.
#' @param xlim 'numeric' vector of length 2.
#' Left and right limits of the spatial grid, data outside these limits is excluded.
#' @param ylim 'numeric' vector of length 2.
#' Lower and upper limits of the spatial grid,
#' data outside these limits is excluded.
#' @param zlim 'numeric' vector of length 2.
#' Minimum and maximum limits of the attribute variable,
#' data outside these limits is excluded.
#' @param ply 'SpatialPolygons', or 'SpatialGridDataFrame'.
#' Cropping polygon
#' @param check_validity 'logical' flag.
#' If true (default), check the validity of polygons.
#' If any of the polygons are invalid, try making them valid by zero-width buffering.
#'
#' @return An object of class 'SpatialPolygonsDataFrame'.
#' The objects \code{data} slot is a data frame, number of rows equal to
#' the number of \code{Polygons} objects and a single column containing attribute values.
#' If \code{level} is true, attribute values are set equal to the midpoint between breakpoints.
#' The status of the polygon as a hole or an island is taken from the ring direction,
#' with clockwise meaning island, and counter-clockwise meaning hole.
#'
#' @note The traditional R graphics model does not draw polygon holes correctly,
#' holes overpaint their containing 'Polygon' object using
#' a user defined background color (white by default).
#' Polygon holes are now rendered correctly using the \code{plot} method for
#' spatial polygons (\code{\link{SpatialPolygons-class}}),
#' see \code{\link{polypath}} for more details.
#' The Trellis graphics model appears to rely on the traditional method so
#' use caution when plotting with \code{\link[sp]{spplot}}.
#'
#' @author J.C. Fisher, U.S. Geological Survey, Idaho Water Science Center
#'
#' @seealso As an alternative, consider using the
#' \code{\link[raster]{rasterToPolygons}} function
#' in the \pkg{raster} package, setting \code{dissolve = TRUE}.
#'
#' @keywords manip
#'
#' @import rgdal
#'
#' @export
#'
#' @examples
#' # Example 1
#' z <- c(1.1, 1.5, 4.2, 4.1, 4.3, 4.7,
#' 1.2, 1.4, 4.8, 4.8, NA, 4.1,
#' 1.7, 4.2, 1.4, 4.8, 4.0, 4.4,
#' 1.1, 1.3, 1.2, 4.8, 1.6, NA,
#' 3.3, 2.9, NA, 4.1, 1.0, 4.0)
#' m <- 5
#' n <- 6
#' x <- rep(0:n, m + 1)
#' y <- rep(0:m, each = n + 1)
#' xc <- c(rep(seq(0.5, n - 0.5, by = 1), m))
#' yc <- rep(rev(seq(0.5, m - 0.5, by = 1)), each = n)
#' grd <- data.frame(z = z, xc = xc, yc = yc, stringsAsFactors = TRUE)
#' sp::coordinates(grd) <- ~ xc + yc
#' sp::gridded(grd) <- TRUE
#' grd <- as(grd, "SpatialGridDataFrame")
#' image(grd, col = gray.colors(30), axes = TRUE)
#' grid(col = "black", lty = 1)
#' points(x = x, y = y, pch = 16)
#' text(cbind(xc, yc), labels = z)
#' text(cbind(x = x + 0.1, y = rev(y + 0.1)),
#' labels = 1:((m + 1) * (n + 1)), cex = 0.6)
#' at <- 1:ceiling(max(z, na.rm = TRUE))
#' plys <- Grid2Polygons(grd, level = TRUE, at = at)
#' cols <- GetColors(length(plys), scheme = "bright", alpha = 0.3)
#' sp::plot(plys, add = TRUE, col = cols)
#' zz <- plys[[1]]
#' legend("top", legend = zz, fill = cols, bty = "n", xpd = TRUE,
#' inset = c(0, -0.1), ncol = length(plys))
#'
#' v1 <- rbind(c( 1.2, 0.5), c(5.8, 1.7), c( 2.5, 5.1), c( 1.2, 0.5))
#' v2 <- rbind(c( 2.5, 2.5), c(3.4, 1.8), c( 3.7, 3.1), c( 2.5, 2.5))
#' v3 <- rbind(c(-0.3, 3.3), c(1.7, 5.1), c(-1.0, 7.0), c(-0.3, 3.3))
#' p1 <- sp::Polygon(v1, hole = FALSE)
#' p2 <- sp::Polygon(v2, hole = TRUE)
#' p3 <- sp::Polygon(v3, hole = FALSE)
#' p <- sp::SpatialPolygons(list(sp::Polygons(list(p1, p2, p3), 1)))
#' plys <- Grid2Polygons(grd, level = TRUE, at = at, ply = p)
#' cols <- GetColors(length(zz), scheme = "bright", alpha = 0.6)
#' cols <- cols[zz %in% plys[[1]]]
#' sp::plot(plys, col = cols, add = TRUE)
#' text(cbind(xc, yc), labels = z)
#'
#' # Example 2
#' data(meuse.grid, package = "sp")
#' sp::coordinates(meuse.grid) <- ~ x + y
#' sp::gridded(meuse.grid) <- TRUE
#' meuse.grid <- as(meuse.grid, "SpatialGridDataFrame")
#' meuse.plys <- Grid2Polygons(meuse.grid, "dist", level = FALSE)
#' op <- par(mfrow = c(1, 2), oma = rep(0, 4), mar = rep(0, 4))
#' sp::plot(meuse.plys, col = heat.colors(length(meuse.plys)))
#' title("level = FALSE", line = -7)
#' meuse.plys.lev <- Grid2Polygons(meuse.grid, "dist", level = TRUE)
#' sp::plot(meuse.plys.lev, col = heat.colors(length(meuse.plys.lev)))
#' title("level = TRUE", line = -7)
#' par(op)
#'
#' # Example 3
#' m <- datasets::volcano
#' m <- m[nrow(m):1, ncol(m):1]
#' x <- seq(from = 2667405, length.out = ncol(m), by = 10)
#' y <- seq(from = 6478705, length.out = nrow(m), by = 10)
#' r <- raster::raster(m, xmn = min(x), xmx = max(x), ymn = min(y),
#' ymx = max(y), crs = "+init=epsg:27200")
#' plys <- Grid2Polygons(r, level = TRUE)
#' cols <- GetColors(length(plys), scheme = "DEM screen")
#' sp::plot(plys, col = cols, border = "#515151")
#'
Grid2Polygons <- function(grd, zcol=1, level=FALSE, at=NULL, cuts=20,
pretty=FALSE, xlim=NULL, ylim=NULL, zlim=NULL,
ply=NULL, check_validity=TRUE) {
# check arguments
stopifnot(inherits(grd, c("BasicRaster",
"SpatialPixelsDataFrame",
"SpatialGridDataFrame")))
checkmate::qassert(zcol, c("S1", "X1[0,)"))
checkmate::assertFlag(level)
checkmate::assertNumeric(at, any.missing=FALSE, min.len=2, null.ok=TRUE)
checkmate::assertCount(cuts, positive=TRUE)
checkmate::assertFlag(pretty)
checkmate::assertNumeric(xlim, len=2, sorted=TRUE, null.ok=TRUE)
checkmate::assertNumeric(ylim, len=2, sorted=TRUE, null.ok=TRUE)
checkmate::assertNumeric(zlim, len=2, sorted=TRUE, null.ok=TRUE)
checkmate::assertClass(ply, "SpatialPolygons", null.ok=TRUE)
checkmate::assertFlag(check_validity)
# convert grid to 'RasterLayer' class
if (!inherits(grd, "RasterLayer"))
grd <- raster::raster(grd, layer=zcol)
# crop grid data using polygon argument
if (!is.null(ply)) {
crs <- raster::crs(grd)
if (is.na(sp::proj4string(ply))) sp::proj4string(ply) <- crs
if (!sp::identicalCRS(grd, ply)) ply <- sp::spTransform(ply, crs)
grd <- raster::crop(grd, ply, snap="out")
}
# crop grid using xlim and ylim
if (is.null(xlim)) xlim <- c(NA, NA)
if (is.null(ylim)) ylim <- c(NA, NA)
if (is.na(xlim[1])) xlim[1] <- raster::xmin(grd)
if (is.na(xlim[2])) xlim[2] <- raster::xmax(grd)
if (is.na(ylim[1])) ylim[1] <- raster::ymin(grd)
if (is.na(ylim[2])) ylim[2] <- raster::ymax(grd)
e <- raster::extent(c(xlim, ylim))
grd <- raster::crop(grd, e, snap="in")
# crop grid using zlim
if (is.null(zlim)) zlim <- c(NA, NA)
if (is.na(zlim[1])) zlim[1] <- min(grd[], na.rm=TRUE)
if (is.na(zlim[2])) zlim[2] <- max(grd[], na.rm=TRUE)
grd[grd[] < zlim[1] | grd[] > zlim[2]] <- NA
grd <- raster::trim(grd)
# determine break points
if (level) {
if (is.null(at)) {
if (pretty)
at <- pretty(zlim, cuts)
else
at <- seq(zlim[1], zlim[2], length.out=cuts)
}
zc <- at[seq_len(length(at) - 1L)] + diff(at) / 2
z <- zc[findInterval(grd[], at, rightmost.closed=TRUE)]
} else {
z <- as.numeric(grd[])
}
# define nodes and elements
m <- dim(grd)[1]
n <- dim(grd)[2]
dx <- raster::xres(grd)
dy <- raster::yres(grd)
xmin <- raster::xmin(grd)
xmax <- raster::xmax(grd)
ymin <- raster::ymin(grd)
ymax <- raster::ymax(grd)
x <- seq(xmin, xmax, by=dx)
y <- seq(ymin, ymax, by=dy)
nnodes <- (m + 1L) * (n + 1L)
nelems <- m * n
nodes <- seq_len(nnodes)
elems <- seq_len(nelems)
coords <- cbind(x=rep(x, m + 1L), y=rep(rev(y), each=n + 1L))
n1 <- unlist(lapply(seq_len(m), function(i) seq(1, n) + (i - 1L) * (n + 1L)))
n2 <- n1 + 1L
n4 <- unlist(lapply(seq_len(m), function(i) seq(1, n) + i * (n + 1L)))
n3 <- n4 + 1L
elem_nodes <- cbind(n1, n2, n3, n4)
# define segments in each element
nsegs <- nelems * 4L
segs <- matrix(data=NA, nrow=nsegs, ncol=4,
dimnames=list(1:nsegs, c("elem", "a", "b", "z")))
segs[, 1] <- rep(seq_len(nelems), each=4)
segs[, 2] <- c(t(elem_nodes))
segs[, 3] <- c(t(elem_nodes[, c(2, 3, 4, 1)]))
segs[, 4] <- rep(z, each=4)
segs <- stats::na.omit(segs)
# identify levels (or unique values)
levs <- sort(unique(stats::na.omit(z)))
# find polygon nodes for each level
poly_nodes <- lapply(levs, function(x) {
FindPolyNodes(segs[segs[, "z"] == x, c("a", "b")])
})
# build 'SpatialPolygons' object
sp_polys <- sp::SpatialPolygons(lapply(seq_along(poly_nodes), function(i) {
sp::Polygons(lapply(poly_nodes[[i]], function(x) {
sp::Polygon(coords[x, ])
}), ID=paste0("X", i))
}), proj4string=raster::crs(grd))
# assign ownership of holes to parent polygons
sp_polys <- rgeos::createSPComment(sp_polys)
# check validity of polygons, if invalid, try making valid by zero-width buffering
if (check_validity) {
if (isFALSE(all(suppressWarnings(rgeos::gIsValid(sp_polys, byid=TRUE)))))
sp_polys <- rgeos::gBuffer(sp_polys, byid=TRUE, width=0)
if (isFALSE(all(rgeos::gIsValid(sp_polys, byid=TRUE))))
stop("invalid polygons found")
}
# convert to 'SpatialPolygonsDataFrame' object, add data frame of levels
d <- data.frame(z=levs, row.names=row.names(sp_polys), stringsAsFactors=TRUE)
sp_polys <- sp::SpatialPolygonsDataFrame(sp_polys, data=d, match.ID=TRUE)
# crop 'SpatialPolygonsDataFrame' object using polygon argument
if (!is.null(ply)) sp_polys <- SetPolygons(sp_polys, ply, "gIntersection")
sp_polys
}
#' Find Polygon Nodes
#'
#' @param s 'matrix'.
#' A 2-column table giving start- and end-node indexes for each segment in a level.
#'
#' @return An object of class 'list'.
#' Vector components giving node indexes for each polygon ring.
#' The status of the polygon as a hole or an island is taken from the ring direction,
#' with clockwise meaning island, and counter-clockwise meaning hole.
#'
#' @author J.C. Fisher, U.S. Geological Survey, Idaho Water Science Center
#'
#' @references \url{https://stackoverflow.com/questions/643995}
#'
#' @keywords internal
#'
#' @useDynLib inlmisc, .registration=TRUE, .fixes="C_"
#'
#' @export
#'
FindPolyNodes <- function(s) {
mode(s) <- "integer"
# remove duplicate segments
id <- paste(pmin.int(s[, 1], s[, 2]), pmax.int(s[, 1], s[, 2]))
s <- s[!id %in% id[duplicated(id)], , drop=FALSE]
# call c program to define polygon rings
out <- matrix(.Call(C_DefinePolygons, s[, 1], s[, 2]), nrow=nrow(s), ncol=2)
# place returned array into list object
poly_nodes <- lapply(unique(out[, 2]), function(x) out[out[, 2] == x, 1])
# close polygon by joining the first point to the last point
poly_nodes <- lapply(poly_nodes, function(x) c(x, x[1]))
poly_nodes
}
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.