CODEWITHSUNDEEP
rr-raster
PaulG
PaulG

Reputation: 53

Speeding up a loop over rasters

I have a big dataset with 30000 rasters. My goal is to extract a mean value using the polygon located within the raster and create a file with extracted rasters values and dates from rasters filenames.

I succeeded in doing this by performing the following loop:

for (i in 1:length(rasters2014)){
a <- raster(rasters2014[i])
ext[i] <- as.vector(extract(a, poligon2, fun=mean, na.rm=TRUE, df=F))
}
output2 = data.frame(ext, filename=filename2014)

The problem is that the presented above loop takes about 2.5h hours to complete the calculation. Does anyone have an idea how I could speed up this process?

Upvotes: 0

Views: 629

Answers (3)

lbusett
lbusett

Reputation: 5932

If your raster are all properly aligned (same ncol, nrow, extent, origin, resolution), you could try identifying the "cell numbers" to be extracted by looking on the first file, then extracting based on those. This could speed-up the processing beacause raster does not need to compute which cells to extract. Something like this: 

rast1 <- raster(rasters2014[1])
cells <- extract(rast1, poligon2, cellnumbers = TRUE, df = TRUE)[,"cells"]
ext <- list()

for (i in 1:length(rasters2014)){
  a <- raster(rasters2014[i])
  ext[[i]] <- as.vector(extract(a, cells, fun=mean, na.rm=TRUE, df=F))
}

Note that I am also using a list to store the results to avoid "growing" a vector, which is usually wasteful.

Alternatively, as suggested by @qdread, you could build a rasterStack using raster::stack(rasters2014, quick = TRUE) and call extract over the stack to avoid the for loop. Don't know which would be faster.

HTH

Upvotes: 2

PaulG
PaulG

Reputation: 53

Thanks for your help! I've tried all of the proposed solutions and the computation time generally the same regardless of the applied method. Therefore, I guess that it is just not possible to significantly speed up the computational process.

Upvotes: 0

Robert Hijmans
Robert Hijmans

Reputation: 47071

If your polygons do not overlap (and in most cases they don't) an alternative route is 

library(raster)
x <- rasterize(poligon2, rasters2014[1])
s <- raster::stack(rasters2014, quick = TRUE)
z <- zonal(s, x, "mean")

PS: Faster is nicer, but I would suggest getting lunch while this runs.

Upvotes: 1

Related Questions