Skip to content

Instantly share code, notes, and snippets.

@JoshOBrien
Last active October 25, 2021 19:22
Show Gist options
  • Save JoshOBrien/7cf19b8b686e6d6230a78a1a9799883b to your computer and use it in GitHub Desktop.
Save JoshOBrien/7cf19b8b686e6d6230a78a1a9799883b to your computer and use it in GitHub Desktop.
gRasterize -- Fast rasterization of spatial objects using gdalUtilities::gdal_rasterize
##' Rasterize \code{Spatial*} objects using gdal_rasterize.
##'
##' For a 1000-by-1000 raster, \code{gRasterize} is more than 6 times
##' faster than \code{raster::rasterize}. For a 2000-by-2000 raster,
##' it is almost 12 times faster (6 seconds vs. 70 seconds on my
##' Windows laptop).
##'
##' I've modeled \code{gRasterize} arguments and behavior on that of
##' \code{\link[raster]{rasterize}}. Like\code{rasterize},it takes a
##' \code{filename=} argument which defaults to \code{""} in which
##' case (unless it's determined internally that
##' \code{!canProcessInMemory(rstr, 3)}) the returned raster is
##' \code{'inMemory'}. Otherwise, if the raster is too large or if a
##' filename is supplied, it's returned \code{fromDisk}.
##'
##' Internally, \code{\link{gdal_rasterize}} by default writes to a
##' file, and only optionally returns an R \code{Raster} object (when
##' its \code{output_Raster = TRUE}); to get the raster
##' \code{'inMemory'}, I use \code{\link[raster]{readAll}} (after a
##' check that it's really OK, memory-wise to do so).
##' @title Fast rasterize for Spatial objects
##' @param SPDF A \code{Spatial*} object to be rasterized.
##' @param r A \code{Raster*} object to be used as the rasterization
##' template.
##' @param field Character. The name of the numeric column in
##' \code{data.frame()} that will be written to the output Raster*
##' object.
##' @param filename Character. Output filename (optional). If none is
##' supplied, the resulting raster will be stored \code{inMemory}
##' (unless it is too large, as determined by a call to
##' \code{canProcessInMemory(y, 3)}
##' @return A \code{RasterLayer} object containing a rasterized
##' version of \code{SPDF}.
##' @importFrom raster levels canProcessInMemory init raster
##' rasterTmpFile ratify readAll shapefile
##' @export
##' @author Joshua O'Brien
##' @examples
##' \donttest{
##' SPDF <- shapefile(system.file("external/lux.shp", package="raster"))
##' ## rr <- raster(extent(SPDF), ncol=100, nrow=100, crs=proj4string(SPDF))
##' llratio <- 1/cos(pi*mean(coordinates(SPDF)[,2])/180)
##' rr <- raster(extent(SPDF),
##' resolution=c(llratio*0.01, 0.01),
##' crs=proj4string(SPDF))
##'
##' ## An example using an integer-valued field
##' rInt <- gRasterize(SPDF, rr, field = "ID_2")
##' plot(rInt, col = RColorBrewer::brewer.pal(name = "Paired", 12))
##' plot(SPDF, lwd = 3, border = "grey30", add = TRUE)
##'
##' ## An example using a character-valued field
##' rFac <- gRasterize(SPDF, rr, field = "NAME_2")
##' rasterVis::levelplot(rFac)
##' }
library(raster)
library(gdalUtilities)
gRasterize <- function (SPDF, r, field, filename = "") {
INMEM <- canProcessInMemory(r, 3) && filename == ""
## character/factor field preparations
ofield <- SPDF[[field]]
FACTOR_FIELD <- inherits(ofield, c("character", "factor"))
if(FACTOR_FIELD) {
SPDF[[field]] <- as.integer(factor(ofield, levels = unique(ofield)))
}
## Write both layers to files ("*.shp" & "*.tif" respectively)
tmpDir <- dirname(rasterTmpFile())
shpDir <- file.path(tmpDir, "shapes")
shpFile <- file.path(shpDir, "SPDF.shp")
if(!dir.exists(shpDir)) {
on.exit(unlink(shpDir, recursive=TRUE), add=TRUE)
dir.create(shpDir)
}
if (filename == "") {
rasFile <- file.path(tmpDir, "r.tif")
} else {
rasFile <- filename
}
shapefile(SPDF, shpFile, overwrite = TRUE)
## We want to burn into a raster with background of NAs
## r[] <- NA
## suppressWarnings(writeRaster(r, filename = rasFile,
## format="GTiff", datatype='FLT4S',
## overwrite = TRUE))
init(r, fun = function(x) NA, filename = rasFile, format = "GTiff",
datatype = "FLT4S", overwrite = TRUE)
## Command-line option processing
## (Will eventually be much more extensive, perhaps in its own function.)
opts <- c("-a", field)
## Rasterize polygons, burning into the raster values from selected field
## b <- gdal_rasterize(shpFile, rasFile, a = field, output_Raster = TRUE)
b <- gdal_rasterize(shpFile, rasFile, a = field)
## r <- raster(b, layer=1)
r <- raster(rasFile)
## Further processing for character or factor fields from SPDF
if(FACTOR_FIELD) {
r <- ratify(r)
RAT <- levels(r)[[1]]
RAT[field] <- unique(ofield)[RAT[[1]]]
levels(r) <- RAT
}
## Clean up intermediate files:
unlink(dirname(shpFile), recursive = TRUE)
if(INMEM) {
r <- readAll(r)
unlink(rasFile)
}
r
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment