I have a rasterBrick/stack named as 'wnt', which has 8 layers. I want to reclassify each layer based on the mean (mn) and standard deviation (sd) of corresponding layer. The classification scheme is like (if the layer's cell values = x) then
x>mn + sd = 1, mn +0.5sd < x < mn+ sd = 2, mn- 0.5sd < x < mn + 0.5sd = 3, mn- sd < x < mn- 0.5sd = 4, x < mn + sd = 5.
The rasterBrick/stack 'wnt' can be found here
I am expecting a final raster stack with 8 reclassified layer.
Following Code I have tried. But not getting the expected result.
# create an empty raster brick to store the reclassified data
wnt_reclass <- brick(nrows=nlayers(wnt), ncols=ncell(wnt),
xmn=extent(wnt)[1], xmx=extent(wnt)[2],
ymn=extent(wnt)[3], ymx=extent(wnt)[4],
crs=projection(wnt))
# loop through each layer of the raster brick, calculate the mean and standard deviation
# for the current layer, and reclassify the data based on the current layer's mean and sd
for (i in 1:nlayers(wnt)) {
# extract the current layer
layer <- wnt[[i]]
# calculate the mean and sd for the current layer
mn <- mean(layer[], na.rm = TRUE)
sd <- sd(layer[], na.rm = TRUE)
# reclassify the data based on the current layer's mean and sd
x <- layer[]
x[x > (mn + sd)] <- 1
x[(mn + 0.5 * sd) < x & x < (mn + sd)] <- 2
x[(mn - 0.5 * sd) < x & x < (mn + 0.5 * sd)] <- 3
x[(mn - sd) < x & x < (mn - 0.5 * sd)] <- 4
x[x < (mn + sd)] <- 5
# set the reclassified data for the current layer in the output raster brick
wnt_reclass[[i]] <- setValues(layer, x)
}
# Set the names of the reclassified layers
names(wnt_reclass) <- paste0("wnt_", 1:nlayers(wnt_reclass))
# Set the extent and projection of the reclassified rasterBrick to match the original 'wnt' rasterBrick
extent(wnt_reclass) <- extent(wnt)
projection(wnt_reclass) <- projection(wnt)
Here is how you can do that with "terra" (the replacement of "raster").
First create some example data (when asking questions, do not include files, except for files that ship with R for use in the examples).
library(terra)
set.seed(1)
r <- rast(nrow=10, ncol=10, nlyr=8, vals=runif(800))
Get the global mean and sd
rmn <- global(r, "mean", na.rm=TRUE)
rsd <- global(r, "sd", na.rm=TRUE)
Write a function that classifies a single layer
f <- function(x, mn, sd) {
classify(x, c(-Inf, mn + sd, mn +0.5*sd, mn-0.5*sd, mn-sd, Inf))
}
Call that function for each layer, using layer specific values for arguments mn and sd
s <- lapply(1:nlyr(r), function(i) f(r[[i]], mn=rmn[i,1], sd=rsd[i,1]))
Combine the list into a SpatRaster and add one to remove the labels, and get a range from 1 to 5 instead of 0 to 4
x <- rast(s) + 1
x
#class : SpatRaster
#dimensions : 10, 10, 8 (nrow, ncol, nlyr)
#resolution : 36, 18 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#source(s) : memory
#names : lyr.1, lyr.2, lyr.3, lyr.4, lyr.5, lyr.6, ...
#min values : 1, 1, 1, 1, 1, 1, ...
#max values : 5, 5, 5, 5, 5, 5, ...