I have a raster stack and I would like to create a new stack by normalizing each layer in the original. I can do it in the inelegant example below. Is there an easier way?
library(raster)
#> Loading required package: sp
set.seed(1)
# Create random data for 3 layers in a stack
s <- stack(raster(matrix(rnorm(100, mean = 10, sd = 1), nrow = 10)),
raster(matrix(rnorm(100, mean = 20, sd = 3), nrow = 10)),
raster(matrix(rnorm(100, mean = 30, sd = 5), nrow = 10)))
# Roll my own functions for normalizing
mu <- function (x) {
mean(extract(x, 1:ncell(x)))
}
sigma <- function (x) {
sd(extract(x, 1:ncell(x)))
}
normalize <- function (x) {
return((x - mu(x))/sigma(x))
}
normalize_stack <- function (x) {
ls <- list()
for (i in 1:dim(x)[3]) {
ls[[i]] <- normalize(subset(x, i))
}
return(stack(ls))
}
# Normalize the stack
normal <- normalize_stack(s)
# Verify that means are essentially 0
for (i in 1:dim(normal)[3]) {
print(mu(subset(normal, i)))
}
#> [1] 7.222304e-16
#> [1] -2.867314e-16
#> [1] -1.519569e-16
# Verify that sds are 1
for (i in 1:dim(normal)[3]) {
print(sigma(subset(normal, i)))
}
#> [1] 1
#> [1] 1
#> [1] 1
Created on 2019-10-16 by the reprex package (v0.3.0)