[R-sig-Geo] raster:::stack(), missing value where TRUE/FALSE needed

Howard, Tim G (DEC) tim.howard at dec.ny.gov
Thu Aug 25 22:44:02 CEST 2016


All,
I'm stymied by this, I expect there is a simple solution. I am trying to stack/brick some very large rasters (26 GB in native R [.gri] format!). But the error does not imply there's a problem with the raster size.  AND, my stack calls work with a crop() of the rasters, as follows :

### first try, has always worked in the past

> library(raster)
> pathToRas <- "G:/test/nativeR"
> pathToBrick <- "G:/test/brick"
> raslist <- list.files(path = pathToRas, pattern = ".grd$")
> gridlist<-as.list(paste(pathToRas,raslist,sep = "/"))
> nm <- substr(raslist,1,nchar(raslist) - 4)
> names(gridlist)<-nm
> gridlist
$canopy_1cell
[1] "G:/test/nativeR/canopy_1cell.grd"

$shrubscrub_10
[1] "G:/test/nativeR/shrubscrub_10.grd"

> envBrick <- brick(stack(gridlist),filename = paste(pathToBrick,"/brick.grd",sep = ""))
Error in if (common.len == 1L) unlist(x, recursive = FALSE) else if (common.len >  : 
  missing value where TRUE/FALSE needed


##### second attempt, load as rasters first

> r1 <- raster(gridlist[[1]])
> r2 <- raster(gridlist[[2]])
> s <- stack(r1, r2)
Error in if (common.len == 1L) unlist(x, recursive = FALSE) else if (common.len >  : 
  missing value where TRUE/FALSE needed


#### third attempt, crop them, write them, then use exactly the same call to create the brick. Works!

> cropExtent <- extent(1950000, 2200000, 2815000, 3050000)
> r1crop <- crop(r1, cropExtent)
> r2crop <- crop(r2, cropExtent)
> writeRaster(r1crop, filename = paste(pathToBrick,"/r1crop", sep=""), format = "raster")
> writeRaster(r2crop, filename = paste(pathToBrick,"/r2crop", sep=""), format = "raster")
> cropgridlist <- as.list(c(paste(pathToBrick,"r1crop", sep="/"), paste(pathToBrick,"r2crop", sep="/")))
> names(cropgridlist) <- c("r1crop","r2crop")
> envBrick <- brick(stack(cropgridlist),filename = paste(pathToBrick,"cropBrick.grd",sep = "/"))
>
> envBrick
class       : RasterBrick 
dimensions  : 6766, 8333, 56381078, 2  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 1950015, 2200005, 2815005, 3017985  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
data source : G:\test \brick\cropBrick.grd 
names       :   r1crop,   r2crop 
min values  :        0,        0 
max values  : 92.55556,  1.00000 

>
#### here's the traceback() if that is useful
10: simplify2array(answer, higher = (simplify == "array"))
9: sapply(1:nl, function(i) {
       r <- x at layers[[i]]
       r at data@names <- value[i]
       r
   })
8: `names<-`(`*tmp*`, value = c("Canopy_1cell", "Shrubscrub_10"))
7: `names<-`(`*tmp*`, value = c("Canopy_1cell", "Shrubscrub_10"))
6: .local(x, ...)
5: stack(.makeRasterList(rlist))
4: stack(.makeRasterList(rlist))
3: .local(x, ...)
2: stack(r1, r2)
1: stack(r1, r2)


> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] raster_2.5-8 sp_1.2-1    

loaded via a namespace (and not attached):
[1] rgdal_1.1-3     parallel_3.3.1  tools_3.3.1     Rcpp_0.12.3     grid_3.3.1      lattice_0.20-33

Thanks in advance for any pointers. 
Tim



More information about the R-sig-Geo mailing list