fix & cleanup PR #286

* `mergeRaster` improvements (#283)
* `splitRaster` improvements (#284)
* moved shared examples to `inst/examples/examples_splitRaster.R`
* rename some arguments to make debugging easier
* don't write raster after doing all cropping; instead crop then write each raster tile in sequence
* cleanup

DESCRIPTION

@@ -67,6 +67,7 @@ Imports:
     raster,
     RColorBrewer,
     shiny,
+    snow,
     stats,
     sp,
     stringi,

NAMESPACE

@@ -411,6 +411,7 @@ importFrom(lubridate,duration)
 importFrom(methods,loadMethod)
 importFrom(parallel,clusterApplyLB)
 importFrom(parallel,clusterEvalQ)
+importFrom(raster,'crs<-')
 importFrom(raster,'extent<-')
 importFrom(raster,Which)
 importFrom(raster,cellFromXY)
@@ -434,11 +435,11 @@ importFrom(raster,nlayers)
 importFrom(raster,nrow)
 importFrom(raster,pointDistance)
 importFrom(raster,raster)
-importFrom(raster,rasterize)
 importFrom(raster,res)
 importFrom(raster,returnCluster)
 importFrom(raster,sampleRegular)
 importFrom(raster,setValues)
+importFrom(raster,writeRaster)
 importFrom(raster,xmax)
 importFrom(raster,xmin)
 importFrom(raster,xres)

NEWS

@@ -2,15 +2,16 @@ Known issues: https://github.com/PredictiveEcology/SpaDES/issues
 
 version 1.1.4.xxxx
 =============
-* spread enhancements: circular spreading, landscape-based functions, allow overlapping events
+* fix bug associated with forthcoming `dplyr` update
+* `spread` enhancements: circular spreading, landscape-based functions, allow overlapping events
 * `rings` argument names closer match to `cir`
 * enhancements to `cir` and addition of new function `rings`
 * performance enhancements in `spread`, `rings`, `cir`
 * `Plot` now clips symbols (generally points) to plotting area. This will allow future "wiping" of plot area.
 * `dev()` returns `dev.cur` invisibly, allowing for finer control of plotting devices
 * add `distanceFromEachPoint`, a multipoint version of `raster::distanceFromPoints`
-* improve `splitRaster` (#276)
-* new function `mergeRaster` (#282)
+* improve `splitRaster` (#276, #284)
+* new function `mergeRaster` (#282, #283)
 * remove `gtools` and `secr` from Imports
 * changed and performance boost to `randomPolygons` (uses `spread` internally now)
 * many more and improved unit tests

R/mergeRaster.R

@@ -26,47 +26,7 @@
 #'
 #' @author Yong Luo and Alex Chubaty
 #'
-#' @examples
-#' require(raster)
-#' # an example with dimensions:
-#' # nrow = 77
-#' # ncol = 101
-#' # nlayers = 3
-#' b <- brick(system.file("external/rlogo.grd", package = "raster"))
-#' r <- b[[1]] # use first layer only
-#' nx <- 3
-#' ny <- 4
-#' y0 <- splitRaster(r, nx, ny) # no buffer
-#' y1 <- splitRaster(r, nx, ny, c(10, 10)) # buffer: 10 pixels along both axes
-#' y2 <- splitRaster(r, nx, ny, c(0.5, 0.5)) # buffer: half the width and length of each tile
-#'
-#' # the original raster:
-#' plot(r) # may require a call to `dev()` if using RStudio
-#'
-#' # the split raster:
-#' layout(mat = matrix(seq_len(nx*ny), ncol = nx, nrow = ny))
-#' plotOrder <- c(4,8,12,3,7,11,2,6,10,1,5,9)
-#' invisible(lapply(y0[plotOrder], plot))
-#'
-#' # can be recombined using `raster::merge`
-#' m0 <- do.call(merge, y0)
-#' all.equal(m0, r) ## TRUE
-#'
-#' m1 <- do.call(merge, y1)
-#' all.equal(m1, r) ## TRUE
-#'
-#' m2 <- do.call(merge, y2)
-#' all.equal(m2, r) ## TRUE
-#'
-#' # or recombine using SpaDES::mergeRaster
-#' n0 <- mergeRaster(y0)
-#' all.equal(n0, r) ## TRUE
-#'
-#' n1 <- mergeRaster(y1)
-#' all.equal(n1, r) ## TRUE
-#'
-#' n2 <- mergeRaster(y2)
-#' all.equal(n2, r) ## TRUE
+#' @example inst/examples/example_splitRaster.R
 #'
 setGeneric("mergeRaster", function(x) {
   standardGeneric("mergeRaster")
@@ -109,5 +69,7 @@ setMethod(
       }
       x[[i]] <- crop(r, extent(xminCut, xmaxCut, yminCut, ymaxCut))
     }
-    return(do.call(raster::merge, x))
+    y <- do.call(raster::merge, x)
+    names(y) <- gsub("_tile[0-9].*$", "", names(x[[1]]))
+    return(y)
 })

R/splitRaster.R

@@ -4,11 +4,16 @@
 #' Split rasters can be recombined using \code{do.call(merge, y)} or \code{mergeRaster(y)},
 #' where \code{y <- splitRaster(x)}.
 #'
-#' @param x   The raster to be split.
+#' This function is parallel-aware, using the same mechanism as used in the \code{raster}
+#' package. Specifically, if you start a cluster using \code{\link{beginCluster}}, then
+#' this function will automatically use that cluster. It is always a good
+#' idea to stop the cluster when finished, using \code{\link{endCluster}}.
 #'
-#' @param nx  The number of tiles to make along the x-axis.
+#' @param r       The raster to be split.
 #'
-#' @param ny  The number of tiles to make along the y-axis.
+#' @param nx      The number of tiles to make along the x-axis.
+#'
+#' @param ny      The number of tiles to make along the y-axis.
 #'
 #' @param buffer  Numeric vector of length 2 giving the size of the buffer along the x and y axes.
 #'                If these values less than or equal to \code{1} are used, this
@@ -17,132 +22,165 @@
 #'                of the number of pixels in each tile (rounded up to an integer value).
 #'                Default is \code{c(0, 0)}, which means no buffer.
 #'
+#' @param path    Character specifying the directory to which the split tiles will be saved.
+#'                If missing, the function creates a subfolder in the current working directory
+#'                based on the raster's name (i.e., using \code{names(x)}).
+#'
 #' @return A list (length \code{nx*ny}) of cropped raster tiles.
 #'
 #' @seealso \code{\link{do.call}}, \code{\link[raster]{merge}}, \code{\link{mergeRaster}}.
 #'
 # igraph exports %>% from magrittr
-#' @importFrom raster crop extent rasterize xmax xmin xres ymax ymin yres
+#' @importFrom parallel clusterApplyLB
+#' @importFrom raster crop 'crs<-' extent getCluster returnCluster writeRaster xmax xmin xres ymax ymin yres
 #' @export
 #' @docType methods
 #' @rdname splitRaster
 #'
 #' @author Alex Chubaty and Yong Luo
 #'
-#' @examples
-#' library(raster)
-#' # an example with dimensions:
-#' # nrow = 77
-#' # ncol = 101
-#' # nlayers = 3
-#' b <- brick(system.file("external/rlogo.grd", package = "raster"))
-#' r <- b[[1]] # use first layer only
-#' nx <- 3
-#' ny <- 4
-#' y0 <- splitRaster(r, nx, ny) # no buffer
-#' y1 <- splitRaster(r, nx, ny, c(10, 10)) # buffer: 10 pixels along both axes
-#' y2 <- splitRaster(r, nx, ny, c(0.5, 0.5)) # buffer: half the width and length of each tile
-#'
-#' # the original raster:
-#' plot(r) # may require a call to `dev()` if using RStudio
-#'
-#' # the split raster:
-#' layout(mat = matrix(seq_len(nx*ny), ncol = nx, nrow = ny))
-#' plotOrder <- c(4,8,12,3,7,11,2,6,10,1,5,9)
-#' invisible(lapply(y0[plotOrder], plot))
-#'
-#' # can be recombined using `raster::merge`
-#' m0 <- do.call(merge, y0)
-#' all.equal(m0, r) ## TRUE
-#'
-#' m1 <- do.call(merge, y1)
-#' all.equal(m1, r) ## TRUE
+#' @example inst/examples/example_splitRaster.R
 #'
-#' m2 <- do.call(merge, y2)
-#' all.equal(m2, r) ## TRUE
-#'
-#' # or recombine using SpaDES::mergeRaster
-#' n0 <- mergeRaster(y0)
-#' all.equal(n0, r) ## TRUE
-#'
-#' n1 <- mergeRaster(y1)
-#' all.equal(n1, r) ## TRUE
-#'
-#' n2 <- mergeRaster(y2)
-#' all.equal(n2, r) ## TRUE
-#'
-setGeneric("splitRaster", function(x, nx, ny, buffer) {
+setGeneric("splitRaster", function(r, nx, ny, buffer, path) {
   standardGeneric("splitRaster")
 })
 
 #' @export
 #' @rdname splitRaster
 setMethod(
   "splitRaster",
-  signature = signature(x = "RasterLayer", nx = "integer",
-                        ny = "integer", buffer = "numeric"),
-  definition = function(x, nx, ny, buffer) {
+  signature = signature(r = "RasterLayer", nx = "integer", ny = "integer",
+                        buffer = "numeric", path = "character"),
+  definition = function(r, nx, ny, buffer, path) {
+    checkPath(path, create = TRUE)
+
     if (length(buffer) > 2) {
       warning("buffer contains more than 2 elements - only the first two will be used.")
       buffer <- buffer[1:2]
     } else if (length(buffer) == 1) {
       buffer <- c(buffer, buffer)
     }
-
     if (buffer[1] < 1) {
-      buffer[1] <- ceiling((buffer[1]*(xmax(x) - xmin(x))/nx)/xres(x))
+      buffer[1] <- ceiling((buffer[1]*(xmax(r) - xmin(r))/nx)/xres(r))
     }
     if (buffer[2] < 1) {
-      buffer[2] <- ceiling((buffer[2]*(ymax(x) - ymin(x))/ny)/yres(x))
+      buffer[2] <- ceiling((buffer[2]*(ymax(r) - ymin(r))/ny)/yres(r))
     }
-    ext <- extent(x)
-    tiles <- vector("list", length = nx*ny)
+
+    ext <- extent(r)
+    extents <- vector("list", length = nx*ny)
     n <- 1L
     for (i in seq_len(nx) - 1L) {
       for (j in seq_len(ny) - 1L) {
-        x0 <- ext@xmin + i*((ext@xmax - ext@xmin) / nx) - buffer[1]*xres(x)
-        x1 <- ext@xmin + (i + 1L)*((ext@xmax - ext@xmin) / nx) + buffer[1]*xres(x)
-        y0 <- ext@ymin + j*((ext@ymax - ext@ymin) / ny) - buffer[2]*yres(x)
-        y1 <- ext@ymin + (j + 1L)*((ext@ymax - ext@ymin) / ny) + buffer[2]*yres(x)
-        tiles[[n]] <- crop(x, extent(x0, x1, y0, y1))
+        x0 <- ext@xmin + i*((ext@xmax - ext@xmin) / nx) - buffer[1]*xres(r)
+        x1 <- ext@xmin + (i + 1L)*((ext@xmax - ext@xmin) / nx) + buffer[1]*xres(r)
+        y0 <- ext@ymin + j*((ext@ymax - ext@ymin) / ny) - buffer[2]*yres(r)
+        y1 <- ext@ymin + (j + 1L)*((ext@ymax - ext@ymin) / ny) + buffer[2]*yres(r)
+        extents[[n]] <- extent(x0, x1, y0, y1)
         n <- n + 1L
       }
     }
+
+    cl <- tryCatch(getCluster(), error = function(e) NULL)
+    on.exit(if (!is.null(cl)) returnCluster())
+
+    croppy <- function(i, e, r, path) {
+      filename <- file.path(path, paste0(names(r), "_tile", i, ".grd"))
+      r_i <- crop(r, e[[i]])
+      crs(r_i) <- crs(r)
+      writeRaster(r_i, filename, overwrite = TRUE)
+      return(raster(filename))
+    }
+
+    tiles <- if (!is.null(cl)) {
+      clusterApplyLB(cl = cl, x = seq_along(extents), fun = croppy, e = extents, r = r, path = path)
+    } else {
+      lapply(X = seq_along(extents), FUN = croppy, e = extents, r = r, path = path)
+    }
+
     return(tiles)
 })
 
 #' @export
 #' @rdname splitRaster
 setMethod(
   "splitRaster",
-  signature = signature(x = "RasterLayer", nx = "numeric", ny = "numeric", buffer = "integer"),
-  definition = function(x, nx, ny, buffer) {
-    return(splitRaster(x, as.integer(nx), as.integer(ny), as.numeric(buffer)))
+  signature = signature(r = "RasterLayer", nx = "numeric", ny = "numeric",
+                        buffer = "integer", path = "character"),
+  definition = function(r, nx, ny, buffer, path) {
+    return(splitRaster(r, as.integer(nx), as.integer(ny), as.numeric(buffer), path))
+})
+
+#' @export
+#' @rdname splitRaster
+setMethod(
+  "splitRaster",
+  signature = signature(r = "RasterLayer", nx = "numeric", ny = "numeric",
+                        buffer = "integer", path = "missing"),
+  definition = function(r, nx, ny, buffer) {
+    return(splitRaster(r, as.integer(nx), as.integer(ny), as.numeric(buffer),
+                       path = file.path(getwd(), names(r))))
+})
+
+#' @export
+#' @rdname splitRaster
+setMethod(
+  "splitRaster",
+  signature = signature(r = "RasterLayer", nx = "numeric", ny = "numeric",
+                        buffer = "numeric", path = "character"),
+  definition = function(r, nx, ny, buffer, path) {
+    return(splitRaster(r, as.integer(nx), as.integer(ny), buffer, path))
+})
+
+#' @export
+#' @rdname splitRaster
+setMethod(
+  "splitRaster",
+  signature = signature(r = "RasterLayer", nx = "numeric", ny = "numeric",
+                        buffer = "numeric", path = "missing"),
+  definition = function(r, nx, ny, buffer) {
+    return(splitRaster(r, as.integer(nx), as.integer(ny), buffer,
+                       path = file.path(getwd(), names(r))))
+})
+
+#' @export
+#' @rdname splitRaster
+setMethod(
+  "splitRaster",
+  signature = signature(r = "RasterLayer", nx = "numeric", ny = "numeric",
+                        buffer = "missing", path = "character"),
+  definition = function(r, nx, ny, path) {
+    return(splitRaster(r, as.integer(nx), as.integer(ny), buffer = c(0, 0), path))
 })
 
 #' @export
 #' @rdname splitRaster
 setMethod(
   "splitRaster",
-  signature = signature(x = "RasterLayer", nx = "numeric", ny = "numeric", buffer = "numeric"),
-  definition = function(x, nx, ny, buffer) {
-    return(splitRaster(x, as.integer(nx), as.integer(ny), buffer))
+  signature = signature(r = "RasterLayer", nx = "numeric", ny = "numeric",
+                        buffer = "missing", path = "missing"),
+  definition = function(r, nx, ny) {
+    return(splitRaster(r, as.integer(nx), as.integer(ny), buffer = c(0, 0),
+                       path = file.path(getwd(), names(r))))
 })
 
 #' @export
 #' @rdname splitRaster
 setMethod(
   "splitRaster",
-  signature = signature(x = "RasterLayer", nx = "numeric", ny = "numeric", buffer = "missing"),
-  definition = function(x, nx, ny) {
-    return(splitRaster(x, as.integer(nx), as.integer(ny), buffer = c(0, 0)))
+  signature = signature(r = "RasterLayer", nx = "integer", ny = "integer",
+                        buffer = "missing", path = "character"),
+  definition = function(r, nx, ny, path) {
+    return(splitRaster(r, nx, ny, buffer = c(0, 0), path))
 })
 
 #' @export
 #' @rdname splitRaster
 setMethod(
   "splitRaster",
-  signature = signature(x = "RasterLayer", nx = "integer", ny = "integer", buffer = "missing"),
-  definition = function(x, nx, ny) {
-    return(splitRaster(x, nx, ny, buffer = c(0, 0)))
+  signature = signature(r = "RasterLayer", nx = "integer", ny = "integer",
+                        buffer = "missing", path = "missing"),
+  definition = function(r, nx, ny) {
+    return(splitRaster(r, nx, ny, buffer = c(0, 0),
+                       path = file.path(getwd(), names(r))))
 })