对于大型 SpatRasters，R 中是否有 terra patch 的替代方案？

Question

总体目标：我正在尝试使用栖息地适宜性地图来运行 Circuitscape。

我已经将数据过滤到大于 0.80 的适宜度，因此我目前有一个包含 0 和 1 单元格的大

SpatRaster

（21373 行和 42746 列）。我需要将补丁 ID 分配给任何触摸单元（8 个方向），以便我可以运行我的连接模型！

我目前的问题，我正在使用terra函数

patches()

，但它目前已经运行了3天，还没有完成，我没有办法（据我所知）跟踪进度。

现在我正在努力保持耐心，让补丁功能运行，以防它需要相当长的时间。

我还尝试了

clump()

包中的

raster

功能，但它也有同样的时间问题。

这是我必须要有耐心的事情还是有其他程序/软件包可以加快这个过程？我有 5 个栅格要查找补丁。

我可以访问 ArcGIS Pro，但如果可能的话更愿意保留在 R 中。

Answer 1

以下示例说明了如何估计在给定硬件上修补该大小的数据所需的时间。请注意，在模拟示例栅格时，我做了某种假设，即 10% 的值不是 NA； patch 函数显然是 O(n^2)； 我预计在我的硬件上运行补丁需要 7 到 8 天。

一种可能合理的前进方式似乎是降低光栅的分辨率。您可以使用此处开发的 mdl 等模型来估计较小尺寸需要多长时间。并找到适当的平衡点。

请注意，运行模拟来获得估计值本身花了我 16 分钟。如果我完成一个较小的较大数据点，它会在接近 1 分钟内完成。你可以在那里看到多项式差异。

target_rows <- 21373
target_cols <- 42746

row_frac <- target_rows/(target_cols+target_rows)

library(tidyverse)
library(terra)
set.seed(42)

# i did up to 13 - reduce to 12 if you arent as patient and want a less accurate estimate :) 
(cols <- 2^(2:13))
(rows <- (2^(2:13)/2))

list_of_rasters <- map2_df(cols,rows,
                           \(co,ro){
                             mult <- ro*co
                             cat("\nMaking size\t", mult , "\trows:\t", ro,"\tcols:\t",co)
                             rb_time <- system.time({
                               
                               r <- rast(nrows=ro, ncols=co, xmin=0)
                               
                               to_pop <- sample.int(n = mult,
                                                    size = mult / 10,
                                                    replace = FALSE)
                               r[to_pop] <- 1L
                               
                             })
                             memsize <-  capture.output(pryr::object_size(wrap(r)))
                             cat("\nMemory use\t",memsize)
                             cat("\nElapsed Seconds making the raster:\t",rb_time[["elapsed"]])
                             tibble(
                               raster=list(r),
                               rows=ro,
                               cols=co,
                               memsize_on_disk=memsize,
                               elapsed = rb_time[["elapsed"]])
                           })

list_of_rasters

patches_list <- map(list_of_rasters$raster,
                    \(r){
                      tm_seconds <- system.time({ pr <- patches(r)})
                      cat("\nElapsed Seconds making the patches:\t",tm_seconds[["elapsed"]])
                      tibble(pr=list(pr),tm_seconds=tm_seconds[["elapsed"]])
                    })
patches_list <- list_rbind(patches_list)

(all_info <- bind_cols(list_of_rasters,
                      patches_list))


(mdl <- lm(tm_seconds~poly(rows*cols,degree=2),data=all_info))
summary(mdl)

want_df <- data.frame(
  rows=target_rows,
  cols=target_cols
)
want_df$tm_seconds <- predict(mdl,newdata=want_df)

(all_info2 <- bind_rows(all_info,
                       want_df))

rows_x_cols <- all_info2$rows*all_info2$cols
plot(rows_x_cols,
     all_info2$tm_seconds)

xtc <- seq(from=0,to=target_cols,length.out=500)
xtr <- seq(from=0,to=target_rows,length.out=500)
xl <- xtc*xtr

lines(x=xl,
      y=predict(mdl,
                newdata=list(
        rows=xtr,
        cols=xtc
      )))

seconds_to_days <- function(x){
  round(x/ (60*60*24),4)
}

# how many days ? 
seconds_to_days(all_info2$tm_seconds)