计算几个不同程度的光栅文件的中位数

Posted 2023-02-25

【中文标题】计算几个不同程度的光栅文件的中位数

【英文标题】：calculate median of several raster files with different extent

【发布时间】：2014-02-19 09:54:08

【问题描述】：

我对 R 很陌生，但我遇到了一个问题，到目前为止我找不到解决方案。 我有一个包含 1000 个光栅文件的文件夹。我必须得到每个单元格所有栅格的中位数。

这些文件包含 NoData 单元（我认为因此它们具有不同的范围） 是否有任何解决方案可以遍历文件夹，将所有文件加在一起以获得中位数？

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

我尝试使用光栅堆栈，但由于范围不同，它不起作用。 感谢您的帮助。

【问题讨论】：

拥有 NoData 和拥有不同的范围是两个不同的东西。范围表示矩形的左、右、上、下边缘的位置，该矩形内的像元可能包含也可能不包含 NoData 值。如果文件确实有不同的范围，您是否只想计算所有文件共享区域的中位数（即所有栅格的交集）？

下一个问题是：栅格的对齐方式是否相同？即光栅边缘是否相互对齐？

你不能：为总范围创建一个栅格，将每个栅格叠加到总范围栅格（相乘），创建一个 rasterStack，然后以中值作为函数运行 calc 吗？跨度>

【参考方案1】：

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

创建一些 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 <- '../***/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，例如 this excellent example。

为此，请调整马赛克签名以及如何调用其参数：

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)

可能远非最佳方法，但希望它有所帮助。

【讨论】：

非常感谢您的帮助！不幸的是，我收到一条错误消息，但我找不到它在哪里：rep(value, times = ncell(x)) 中的错误：'times' 参数无效另外：警告消息：在 setValues(x, rep(value, times = ncell(x))) : NAs 由 .local(x, i, j, ..., value) 中的强制错误引入：无法替换此栅格上的值（任何想法都太大了？再次感谢

你是在最开始# Write them out导出r2和r3吗？我刚写了第一个的代码...

@user3327377 你能提供一些你的栅格吗？我不明白你是在你的数据上运行还是在复制上面的代码。

我正在尝试在我自己的数据上运行它。

@user3327377 你能分享一些来自 Dropbox 或 googledrive 的栅格吗？