Search code examples
rstackmeanraster

What is the difference between mean(x) and calc(x, fun=mean) for RasterStack?


I have around 80 raster stacks with 24 layers each. I'm trying to calculate the mean of all layers for each stack (e.g. mean pixel-wise in 24 layers).

I found two approaches to accomplish this but I get slightly different results in the min-max values. Here is an example:

library(raster)
set.seed(9999)

##create example raster list
rasterls=list()
for (i in 1:24){
  r <- raster(nrow=(255*5), ncol=(487*5), vals=runif(255*5*487*5))
  rasterls[[i]] <- r
}
##create raster stack
stackRas <- stack(rasterls)

##calculate mean using mean()
mean(stackRas, na.omit=T)

'''class      : RasterLayer 
dimensions : 1275, 2435, 3104625  (nrow, ncol, ncell)
resolution : 0.1478439, 0.1411765  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 0.2453108, 0.8256463  (min, max)
'''

##calculate mean using calc(fun=mean)
calc(stackRas, fun = mean, na.omit=T)

'''
class      : RasterLayer 
dimensions : 1275, 2435, 3104625  (nrow, ncol, ncell)
resolution : 0.1478439, 0.1411765  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : C:/Users/--/AppData/Local/Temp/RtmpoVVABX/raster/r_tmp_2021-04-28_120400_4844_87419.grd 
names      : layer
values     : 0.2138654, 0.8183816  (min, max)
'''

What is the difference between both approaches? Why the different min-max? Am I missing something? Which one should I use? I'd appreciate any insight. Thanks


Solution

  • There is no na.omit argument for mean. The correct argument should be na.rm. When you passed a value for na.omit, it was ignored and the two methods differ in their default handling of missing values. If we use the correct argument, the results are the same.

    mean(stackRas, na.rm=T)
    # class      : RasterLayer 
    # dimensions : 1275, 2435, 3104625  (nrow, ncol, ncell)
    # resolution : 0.1478439, 0.1411765  (x, y)
    # extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    # crs        : +proj=longlat +datum=WGS84 +no_defs 
    # source     : memory
    # names      : layer 
    # values     : 0.2113316, 0.7704163  (min, max)
    
    calc(stackRas, mean, na.rm=T)
    # class      : RasterLayer 
    # dimensions : 1275, 2435, 3104625  (nrow, ncol, ncell)
    # resolution : 0.1478439, 0.1411765  (x, y)
    # extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    # crs        : +proj=longlat +datum=WGS84 +no_defs 
    # source     : /tmp/RtmpRMo1uo/raster/r_tmp_2021-04-28_124738_12579_61131.grd 
    # names      : layer 
    # values     : 0.2113316, 0.7704163  (min, max)