0

I have 2 raster stacks. I want to do some math between pairs of raster layers from each of the stacks, producing a 3rd raster stack of the same number of layers, ie;

r1<- raster stack 1 # 10 raster layers
r2<- raster stack 2 # 10 raster layers

r3<- sqrt(r1^2 + r2^2) # 10 raster layers

Is this the equivalent of the loop form (for illustrative purposes);

for (i in 1:10) {
r <- sqrt(r1[[i]]^2 + r2[[i]]^2)
r3 <-stack(r3,r)
}

Or is there a more efficient function or apply solution? Thanks.

Jaap
  • 81,064
  • 34
  • 182
  • 193
user2175481
  • 147
  • 1
  • 8
  • 1
    Why not just `r3 <- sqrt((r1 ** 2) + (r2 ** 2))`? – Mislav Jan 29 '18 at 21:17
  • Yes, I think that's what I was asking. So this is the same as the loop I showed taking each pair of raster layers? I wondered if if had to unstack the raster stacks first. – user2175481 Jan 29 '18 at 21:24

1 Answers1

0

Example data

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster")) 

You can indeed do

x <- sqrt(s^2 + s^2)

An alternative approach would be

y <- overlay(s, s, fun=function(x, y) { sqrt(x^2 + y^2) } )

The result is the same:

all(values(x-y) == 0)
#[1] TRUE

A difference is that if you use overlay you can provide a filename to write the result to disk

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63