Iterate through multiple raster stacks

[Iara Lacher]

Hi All,

Does anyone know how to iterate through multiple raster stacks? I need to apply the function to the raster stack itself, not each layer within a stack. I've got ~ 1800 iterations, which is way to much to do one at a time.

Can I use lapply, or a loop, or write a function to do this? Everything I've tried doesn't seem to work. I have looked everywhere online for help.

Thanks a million.

--

Dr. Iara Lacher

Fulbright Scholar - Universidade Federal de Goiás, Brasil
PhD - UC Davis Ecology
MESM - Bren School for Environmental Science and Management

*~* Be a part of the solution ~ Please print only when necessary *~*

[Eric Holmes]

Hi Lara,

I'll preface with fact that I don't use raster stacks often, but I don't see an problem with applying a function to a list of stacks using a loop or lapply.  I have included an example (below) using both methods applying a polygon extract using two stacks of ascii rasters.  Both seem to work just fine although capturing the outputs would be slightly different between methods.  Maybe your function needs to be tweaked.

For large sets of rasters I often try to cut the processing time down by running individual raster operations in parallel with the foreach and doParallel packages, but it sounds like that may not be an option for you.  Hope this helps!

Eric

library(raster)

## Loading required package: sp

library(maptools)

## Checking rgeos availability: TRUE

##List Rasters
rastlist<-list.files("C:/Users/Eric/Desktop/RCAT/Data/PRISM_rasters", full.names=T)
##Create stacks and set projection
s1 <- stack(rastlist[1:10])

## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
## Path to GDAL shared files: C:/Users/Eric/Documents/R/win-library/3.1/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/Eric/Documents/R/win-library/3.1/rgdal/proj

s2 <- stack(rastlist[11:20])
projection(s1) <- "+proj=longlat +datum=NAD83"
projection(s2) <- "+proj=longlat +datum=NAD83"

##read in a zone shapefile for extraction
z <- readShapePoly("C:/Users/Eric/Documents/USFWS/R_Tools/PRISM_Extractor/zones/TJ_Slough_NWR_Boundary.shp", proj4string=CRS("+proj=longlat +datum=NAD83"))

##Extract via loop
for(i in c(s1,s2)){
print(extract(i, z, na.rm=T, fun=mean, weights = TRUE))
}

##      us_ppt_1895.01 us_ppt_1895.02 us_ppt_1895.03 us_ppt_1895.04
## [1,]          18531           2131           3807          921.2
##      us_ppt_1895.05 us_ppt_1895.06 us_ppt_1895.07 us_ppt_1895.08
## [1,]            400            4.8              0              0
##      us_ppt_1895.09 us_ppt_1895.10
## [1,]          446.4          263.2
##      us_ppt_1895.11 us_ppt_1895.12 us_ppt_1895.14 us_ppt_1896.01
## [1,]           3365           1024          30959           4783
##      us_ppt_1896.02 us_ppt_1896.03 us_ppt_1896.04 us_ppt_1896.05
## [1,]          901.6           1328          947.6            246
##      us_ppt_1896.06 us_ppt_1896.07
## [1,]            4.8           11.6

##Extract via lapply
l1 <- list(s1,s2)
lapply(l1, extract, z, na.rm=T, fun=mean, weights = TRUE)

## [[1]]
##      us_ppt_1895.01 us_ppt_1895.02 us_ppt_1895.03 us_ppt_1895.04
## [1,]          18531           2131           3807          921.2
##      us_ppt_1895.05 us_ppt_1895.06 us_ppt_1895.07 us_ppt_1895.08
## [1,]            400            4.8              0              0
##      us_ppt_1895.09 us_ppt_1895.10
## [1,]          446.4          263.2
## 
## [[2]]
##      us_ppt_1895.11 us_ppt_1895.12 us_ppt_1895.14 us_ppt_1896.01
## [1,]           3365           1024          30959           4783
##      us_ppt_1896.02 us_ppt_1896.03 us_ppt_1896.04 us_ppt_1896.05
## [1,]          901.6           1328          947.6            246
##      us_ppt_1896.06 us_ppt_1896.07
## [1,]            4.8           11.

--
Check out our R resources at http://www.noamross.net/davis-r-users-group.html
---
You received this message because you are subscribed to the Google Groups "Davis R Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to davis-rug+...@googlegroups.com.
Visit this group at http://groups.google.com/group/davis-rug.
For more options, visit https://groups.google.com/d/optout.

--

Eric Holmes

Research Ecologist

Center for Watershed Sciences

University of California, Davis

ejho...@ucdavis.edu
(530)305-1688