计算不同程度的几个栅格文件的中位数 [英] calculate median of several raster files with different extent

问题描述

我对R还是很陌生，但有一个问题，到目前为止我还没有找到解决方案. 我有1000个栅格文件的文件夹.我必须获取每个像元的所有栅格的中位数.

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.

文件包含NoData单元(我认为因此它们具有不同的范围) 有什么解决方案可以遍历文件夹，将所有文件加在一起得到中位数?

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.

推荐答案

我将尝试通过mosaic()处理具有不同范围和来源但具有相同分辨率的图像来解决此问题.

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

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，例如为此，请调整镶嵌签名以及如何使用以下命令调用其参数:

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.

最后，函数:

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.

这篇关于计算不同程度的几个栅格文件的中位数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！