вычислить медиану нескольких растровых файлов разной степени

Я попытаюсь подойти к этому с помощью mosaic() изображений с разным размером и происхождением, но с одинаковым разрешением.

Создайте несколько объектов rasterLayer и экспортируйте их (чтобы прочитать последний)

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)

Чтение и мозаика с функцией

Поскольку raster::mosaic не принимает rasterStack/rasterBrick или списки rasterLayers, лучшим подходом является использование do.call, например отличный пример.

Для этого настройте сигнатуру мозаики и способ вызова ее аргументов с помощью:

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)
          })

Давайте сохраним низкую терпимость здесь, чтобы оценить любое неправильное поведение нашей функции.

Наконец, функция:

Функция мозаики

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)

Возможно, далеко не лучший подход, но надеюсь, что это поможет.