[R-sig-Geo] problem with overlay applied to rasterStack

[Melanie Weynants]

Hi,
I am trying to apply a function with overlay on a raster stack. It works OK
with a small stack of synthetic data, but crashes when I apply it to a
larger stack.
Here is my script:

# create test data with structure similar to real data
# 100 x 100 x 7 array of integers between 1 and 4
a <- array(round(runif(70000,1,4)),dim=c(100,100,7))
# add no data/desert/permafrost/water
smpl0 <- sample.int(100,100, replace=TRUE)
a[smpl0[1:50],smpl0[51:100],] <- sample(c(0,5:7), 100, replace=TRUE)
# make brick
b <- brick(a)
# or stack
s <-
stack(raster(a[,,1]),raster(a[,,2]),raster(a[,,3]),raster(a[,,4]),raster(a[,,5]),raster(a[,,6]),raster(a[,,7]))

f <-function(sq1,sq2,sq3,sq4,sq5,sq6,sq7){
   indna <- ! sq1 %in% 1:4 | is.na(sq1)
   sr <- rep(NA, length(sq1))
   sr[indna] <- sq1[indna]
   srlow <- function(sq1,sq3,sq4,sq5,sq6,sq7){
     cb <- cbind(sq4,sq5,sq6,sq7)
     s <- apply(cb, FUN=sort, MARGIN=1)
     means <- rowMeans(t(s)[,-1]) # mean of other 3
     return(sq1*sq3*s[1,]*means)}
   tr <- function(sq){
     sqt <- rep(NA,length(sq))
     sqt[sq==1] <- 0.9
     sqt[sq==2] <- 0.7
     sqt[sq==3] <- 0.5
     sqt[sq==4] <- 0.3
     return(sqt)
   }
   sr[!indna] <- srlow(tr(sq1[!indna]), tr(sq3[!indna]), tr(sq4[!indna]),
tr(sq5[!indna]), tr(sq6[!indna]), tr(sq7[!indna]))

   return(sr)}

# on test data
# works on brick
tmp <- overlay(b, fun=f, filename="tmp.tif",format="GTiff")
# and on stack
tmp <- overlay(s,fun=f)
# until here everything works fine

# but when I want to apply to my data:
download.file("
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/Soil_Quality/sq1.asc",
destfile="sq1.asc", method = "internal", mode="w")
download.file("
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/Soil_Quality/sq2.asc",
destfile="sq2.asc", method = "internal", mode="w")
download.file("
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/Soil_Quality/sq3.asc",
destfile="sq3.asc", method = "internal", mode="w")
download.file("
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/Soil_Quality/sq4.asc",
destfile="sq4.asc", method = "internal", mode="w")
download.file("
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/Soil_Quality/sq5.asc",
destfile="sq5.asc", method = "internal", mode="w")
download.file("
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/Soil_Quality/sq6.asc",
destfile="sq6.asc", method = "internal", mode="w")
download.file("
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/Soil_Quality/sq7.asc",
destfile="sq7.asc", method = "internal", mode="w")

SQ <-
stack("sq1.asc","sq2.asc","sq3.asc","sq4.asc","sq5.asc","sq6.asc","sq7.asc")
projection(SQ) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0"

#class       : RasterStack
#dimensions  : 2160, 4320, 9331200, 7  (nrow, ncol, ncell, nlayers)
#resolution  : 0.08333333, 0.08333333  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0
#names       :         sq1,         sq2,         sq3,         sq4,
sq5,         sq6,         sq7
#min values  :           0, -2147483648, -2147483648, -2147483648,
  0, -2147483648, -2147483648
#max values  :           7,  2147483647,  2147483647,  2147483647,
  7,  2147483647,  2147483647

# it crashes...
SRlow <- overlay(SQ, fun=f, filename = "SRlow.tif")
# Error in .overlayList(x, fun = fun, filename = filename, ...) :
  cannot use this formula, probably because it is not vectorized
# But the function is vectorized and works well with a smaller rasterLayers
object.

# when debugging, it shows that overlay is not running function f on chunks
in all layers but rather it takes one single cell values of each layer and
feeds them as argument sq1 to f, without assigning any values to the other
arguments.

# The following other overlay does not produce any error:
majorSoilConstraint <- function(sq1,sq2,sq3,sq4,sq5,sq6,sq7){
  indna <- ! sq1 %in% 1:4 | is.na(sq1)
  sr <- rep(NA, length(sq1))
  sr[indna] <- sq1[indna]
  constr <- cbind(sq1,sq2,sq3,sq4,sq5,sq6,sq7)
  indnc <- apply(constr[!indna,], FUN=function(x){max(x) == 1}, MARGIN=1) #
no constraints
  sc <- rep(NA, length(sq1))
  sc[indna] <- -sq1[indna]
  sc[!indna][indnc] <- 8
  sc[!indna][!indnc] <- apply(constr[!indna,][!indnc,],
FUN=function(x){which(x==max(x))[1]}, MARGIN=1)
  return(sc)
}

msc <- overlay(SQ, fun=majorSoilConstraint, filename = "msc.tif",
projection="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0", overwrite=TRUE)

# R version 3.2.1 (2015-06-18)
# packageVersion("raster")
#[1] ‘2.4.20’

	[[alternative HTML version deleted]]