calculate median of several raster files with different extent

Question

I'm quite new to R and I have a problem on which I couldn't find a solution so far. I have a folder of 1000 raster files. I have to get the median of all rasters for each cell.

The files contain NoData Cells (I think therefore they have different extents) Is there any solution to loop through the folder, adding together all files an getting the median?

Error in rep(value, times = ncell(x)) : invalid 'times' argument
In addition: Warning message:
In setValues(x, rep(value, times = ncell(x))) : NAs introduced by coercion
Error in .local(x, i, j, ..., value) : 
  cannot replace values on this raster (it is too large

I tried with raster stack, but it doesn't work because of the different extents.
Thanks for your help.

Answer 1

I'll try to approach this by mosaic()'ing images with different extents and origins but same resolution.

Create a few rasterLayer objects and export them (to read latter)

library('raster')
library('rgdal')

e1 <- extent(0,10,0,10)
r1 <- raster(e1)
res(r1) <- 0.5
r1[] <- runif(400, min = 0, max = 1)
#plot(r1)

e2 <- extent(5,15,5,15)
r2 <- raster(e2)
res(r2) <- 0.5
r2[] <- rnorm(400, 5, 1)
#plot(r2)

e3 <- extent(18,40,18,40)
r3 <- raster(e3)
res(r3) <- 0.5
r3[] <- rnorm(1936, 12, 1)
#plot(r3)

# Write them out
wdata <- '../Stackoverflow/21876858' # your local folder
writeRaster(r1, file.path(wdata, 'r1.tif'),
            overwrite = TRUE)
writeRaster(r2,file.path(wdata, 'r2.tif'),
            overwrite = TRUE)
writeRaster(r3,file.path(wdata, 'r3.tif'),
            overwrite = TRUE)

Read and Mosaic'ing with function

Since raster::mosaic do not accept rasterStack/rasterBrick or lists of rasterLayers, the best approach is to use do.call, like this excellent example.

To do so, adjust mosaic signature and how to call its arguments with:

setMethod('mosaic', signature(x='list', y='missing'), 
          function(x, y, fun, tolerance=0.05, filename=""){
            stopifnot(missing(y))
            args <- x
            if (!missing(fun)) args$fun <- fun
            if (!missing(tolerance)) args$tolerance<- tolerance
            if (!missing(filename)) args$filename<- filename
            do.call(mosaic, args)
          })

Let's keep tolerance low here to evaluate any misbehavior of our function.

Finally, the function:

Mosaic function

f.Mosaic <- function(x=x, func = median){
  files <- list.files(file.path(wdata), all.files = F)
  # List  TIF files at wdata folder
  ltif <- grep(".tif$", files, ignore.case = TRUE, value = TRUE) 
  #lext <- list()
  #1rt <- raster(file.path(wdata, i),
  #            package = "raster", varname = fname, dataType = 'FLT4S')
  # Give an extent area here (you can read it from your first tif or define manually)
  uext <- extent(c(0, 100, 0, 100)) 
  # Get Total Extent Area
  stkl <- list()
  for(i in 1:length(ltif)){
    x <- raster(file.path(wdata, ltif[i]),
                package = "raster", varname = fname, dataType = 'FLT4S')
    xext <- extent(x)
    uext <- union(uext, xext)
    stkl[[i]] <- x
  }
  # Global Area empty rasterLayer
  rt <- raster(uext)
  res(rt) <- 0.5
  rt[] <- NA
  # Merge each rasterLayer to Global Extent area
  stck <- list()
  for(i in 1:length(stkl)){
    merged.r <- merge(stkl[[i]], rt,  tolerance = 1e+6)
    #merged.r <- reclassify(merged.r, matrix(c(NA, 0), nrow = 1))
    stck[[i]] <- merged.r
  }
  # Mosaic with Median
  mosaic.r <- raster::mosaic(stck, fun = func) # using median
  mosaic.r
}
# Run the function with func = median
mosaiced <- f.Mosaic(x, func = median)
# Plot it
plot(mosaiced)

Possibly far from the best approach but hope it helps.