来自大型栅格层的随机采样答案

【问题标题】：Random sampling from large rasterlayer来自大型栅格层的随机采样
【发布时间】：2015-02-13 11:52:51
【问题描述】：

我有一个大的光栅层，整数范围从 0 到 44。

class       : RasterLayer
dimensions  : 29800, 34470, 1027206000  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 331300, 676000, 5681995, 5979995  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
data source : /home/mkoehler/stk_rast_whz
names       : stk_rast_whz
values      : 0, 44  (min, max)

我想对每层 5000 个点进行分层抽样。我收到以下错误：

POINTS<-sampleStratified(b, size=5000, na.rm=T, xy=F)
(Error in ys[[i]] <- y : attempt to select less than one element)

这是重现问题的代码（即使只选择 1 每层项目）：

 set.seed(10)
 r <- raster(ncol=5000, nrow=5000)
 names(r) <- 'stratum'
 r[] <- round((runif(ncell(r)))*44)

 sampleStratified(r, size=1,xy=T)

Error in ys[[i]] <- y : attempt to select less than one element

尝试使用更少的层并更改“大小”的设置或 “exp”没有效果。 R 版本：[64 位] C:\Program Files\R\R-3.1.1

有什么想法吗？

提前致谢！

【问题讨论】：

您的示例代码对我来说运行良好（R 3.1.0，32 位）
奇怪！这是否取决于我使用的 R 版本/库版本？
可能是。检查raster 包的版本。我正在使用 2.2-31。
@koekenbakker 我的光栅版本是 raster_2.3-12 我使用相同版本的同事收到相同的错误消息。我会用旧版本试一试。
是的，这是个好主意。您也可以尝试旧版本的代码here。也许您可以通过逐步运行代码并检查每个步骤的输出来进行调试。

标签： r raster

【解决方案1】：

这似乎是一个错误（如raster 2.3-12），当 (1) 栅格包含值为 0 的像元，以及 (2) 栅格无法在内存中处理 (即canProcessInMemory(r) 是FALSE)。

该函数循环遍历freq(r) 生成的唯一单元格值，然后依次按这些值中的每一个索引一个列表。如果这些值之一为零，则将触发错误，因为第 0 个元素不存在。例如：

list()[[0]]
# Error in list()[[0]] : attempt to select less than one element]

您会注意到，如果您用 r[] <- sample(44, ncell(r), replace=TRUE) 等填充 r，则不会发生错误，因为它不会有任何零。

当光栅可以在内存中被处理时，函数循环freq(r)的行号，因此后续的列表索引是明智的。

我已联系维护者报告此错误。

同时，作为临时修复，您可以使用以下内容来制作函数的更正副本（在当前 R 会话中仍然可用）。

sampleStratified2 <- 
  eval(parse(text=sub('sr\\[, 2\\] == i', 'sr[, 2] == f[i, 1]',
                      sub('i in f\\[, 1\\]', 'i in seq_len(nrow(f))',
                          deparse(getMethod(sampleStratified, 
                                            signature='RasterLayer')@.Data))
  )))

sampleStratified2(r, size=1, xy=TRUE)

【讨论】：

谢谢，这已在光栅版本 2.3-20 中得到修复（目前仅可从 R-Forge 获得，即将在 CRAN 上使用）。
@RobertH 太好了 - 感谢您这么快解决这个问题。