2
\$\begingroup\$

I have several hundred "A" rasters, numbered from 001_A.tif to 100_A.tif.

I also have several hundred "B" rasters, numbered from 001_B.tif to 100_B.tif.

I want to multiply them like this:

001_A * 001_B = 001_result
002_A * 002_B = 002_result
003_A * 003_B = 003_result
.
.
.

I have a loop in R which works fine, but it is relatively slow:

library(raster)

A_Files <- list.files(pattern="A.tif") # files A
All_A <- unique(A_Files)

B_Files <- list.files(pattern="B") # files B
All_B <- unique(B_Files)

for (i in 1:100) {
  A.i <- All_A[i]
  B.i <- All_B[i]

  Raster_A.i <- raster(A.i) # create rasters from files A
  Raster_B.i <- raster(B.i) # create rasters from files B

  Result.i <- (Raster_A.i * Raster_B.i) # multiplying A by B

  writeRaster(Result.i, paste(A.i, "_result.tif", sep = ""), datatype='FLT4S',
 overwrite=TRUE) # write results in .tif format with the name "001_A_result.tif"
    } 

Is there a way to speed it up?

\$\endgroup\$
1

2 Answers 2

1
\$\begingroup\$

Rasters

Rasters are highly memory intensive in R's already fragile memory space. However, one does not necessarily need to read a raster into memory each time. In fact, it is highly frowned upon to read in a raster in that way as it just eats up so much memory and is computationally intensive.

stack()'s beauty

It's much better to make a stack of rasters as it does not load them into memory, but opens their file location.

library(raster)

A_Files <- list.files(pattern="A.tif") # files A
All_A <- stack(unique(A_Files))

B_Files <- list.files(pattern="B") # files B
All_B <- stack(unique(B_Files))

Notice, all of the elements are now in a "vectorized" form.

overlay() to calculate

From here, we can use overlay() to perform a calculation. Note, since the number of layers is equivalent between (ALL_A) and (ALL_B) we do not need to worry about an element being repeated.

x <- overlay(All_A, All_B, fun=function(x,y) x*y )

Export

The only downside to this approach is you will have to go through a loop to export the answers:

out <- unstack(x)
outputnames <- paste0(seq_along(out), "_result.tif")
for(i in seq_along(out)){ writeRaster(out[[i]], file=outputnames[i]) }

References

For multiplying stacks see:

For more information on stacks (see bricks as well)

\$\endgroup\$
1
  • \$\begingroup\$ Thank you for your answer. But the overlay part is not working. Or maybe it is working, but it takes so long that I interrupted it manually. I chcecked it by the system.time function. Doesn't it multiply each raster A with each raster B, instead of multiplying only two matching rasters? \$\endgroup\$ Feb 22, 2016 at 14:26
1
\$\begingroup\$

Finally I saved some time using the lapply or the map function.

Here is the code with lapply:

require(raster) # raster package

setwd("C:/Rasters") # set working directory

A_Files <- list.files(pattern="A.tif") # files A
All_A <- stack(unique(A_Files)) # stacking As

B_Files <- list.files(pattern="B") # files B
All_B <- stack(unique(B_Files)) # stacking Bs

As <- lapply(All_A, raster) # rasters from A

Bs <- lapply(All_B, raster) # rasters from B

# multiplying A*B:
Multiply_A_B <- lapply(seq_along(As), FUN = function(i) As[[i]] * Bs[[i]])

# saving results as "_result.tif"
lapply(seq_along(Multiply_A_B), function(x) {
    writeRaster(Multiply_A_B[[x]],
    paste(All_A[x], "_result.tif", sep = ""),
      datatype = 'FLT4S', overwrite = TRUE)
})

Code with the map function:

# multiplying A*B:
Multiply <- Map("*", As, Bs)
# the rest is the same as with the lapply example above

The operation took 210 seconds with for loop, while with the lapply it took 190 seconds and with map it took 185 seconds.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.