select 栅格数据有效范围的最快方法

Question

使用 R，我需要以最快的方式 select 给定栅格（来自包 raster）的有效范围。我试过这个：

library(raster)
library(microbenchmark)
library(ggplot2)
library(compiler)

r <- raster(ncol=100, nrow=100)
r[] <- runif(ncell(r))

#Let's see if precompiling helps speed...
f <- function(x, min, max) reclassify(x, c(-Inf, min, NA, max, Inf, NA))
g <- cmpfun(f)

#Benchmark!
compare <- microbenchmark(
    calc(r, fun=function(x){ x[x < 0.2] <- NA; x[x > 0.8] <- NA; return(x)}), 
    reclassify(r, c(-Inf, 0.2, NA, 0.8, Inf, NA)),
    g(r, 0.2, 0.8),
    times=100)
autoplot(compare) #Reclassify is much faster, precompiling doesn't help much.

#Check they are the same...
identical(
          calc(r, fun=function(x){ x[x < 0.2] <- NA; x[x > 0.8] <- NA; return(x)}),
          reclassify(r, c(-Inf, 0.2, NA, 0.8, Inf, NA))
) #TRUE
identical(
          reclassify(r, c(-Inf, 0.2, NA, 0.8, Inf, NA)),
          g(r, 0.2, 0.8),
) #TRUE

重新分类方法要快得多，但我确信它可以更快。我该怎么做？

Answer 1

还有一种方法：

h <- function(r, min, max) {
  rr <- r[]
  rr[rr < min | rr > max] <- NA
  r[] <- rr
  r
}

i <- cmpfun(h)

identical(
  i(r, 0.2, 0.8),
  g(r, 0.2, 0.8)
)



#Benchmark!
compare <- microbenchmark(
  calc(r, fun=function(x){ x[x < 0.2] <- NA; x[x > 0.8] <- NA; return(x)}), 
  reclassify(r, c(-Inf, 0.2, NA, 0.8, Inf, NA)),
  g(r, 0.2, 0.8),
  h(r, 0.2, 0.8),
  i(r, 0.2, 0.8),
  times=100)
autoplot(compare)

在这种情况下编译没有太大帮助。

您甚至可以通过直接使用 @ 访问光栅对象的插槽来进一步加快速度（尽管通常不鼓励这样做）。

j <- function(r, min, max) {
  v <- r@data@values
  v[v < min | v > max] <- NA
  r@data@values <- v
  r
}

k <- cmpfun(j)

identical(
  j(r, 0.2, 0.8)[],
  g(r, 0.2, 0.8)[]
)

Answer 2

光栅包有一个函数：clamp。它比 g 快，但比 h 和 i 慢，因为它内置了一些开销（安全）。

compare <- microbenchmark(
  h(r, 0.2, 0.8),
  i(r, 0.2, 0.8),
  clamp(r, 0.2, 0.8),
  g(r, 0.2, 0.8),
  times=100)
autoplot(compare)

Answer 3

虽然这个问题的可接受答案对于示例栅格是正确的，但重要的是要注意最快的安全函数在很大程度上取决于栅格大小：函数 h 和 i 提供@rengis 只有在使用相对较小的栅格（以及相对简单的重新分类）时才会更快。只需将 OP 示例中光栅 r 的大小增加十倍，即可使 reclassify 更快：

# Code from OP @AF7
library(raster)
library(microbenchmark)
library(ggplot2)
library(compiler)

#Let's see if precompiling helps speed...
f <- function(x, min, max) reclassify(x, c(-Inf, min, NA, max, Inf, NA))
g <- cmpfun(f)

# Funcions from @rengis
h <- function(r, min, max) {
  rr <- r[]
  rr[rr < min | rr > max] <- NA
  r[] <- rr
  r
}

i <- cmpfun(h)

# Benchmark with larger raster (100k cells, vs 10k originally)
r <- raster(ncol = 1000, nrow = 100)
r[] <- runif(ncell(r))

compare <- microbenchmark(
  calc(r, fun=function(x){ x[x < 0.2] <- NA; x[x > 0.8] <- NA; return(x)}), 
  reclassify(r, c(-Inf, 0.2, NA, 0.8, Inf, NA)),
  g(r, 0.2, 0.8),
  h(r, 0.2, 0.8),
  i(r, 0.2, 0.8),
  times=100)
autoplot(compare)

reclassify 变得更快的确切点取决于栅格中的像元数量和重新分类的复杂性，但在这种情况下，交叉点大约为 50,000 个像元（见下文）。

随着栅格变得更大（或计算更复杂），另一种加速重分类的方法是使用多线程，例如使用 snow 包：

# Reclassify, using clusterR to split into two threads
library(snow)
tryCatch({
      beginCluster(n = 2)
      clusterR(r, reclassify, args = list(rcl = c(-Inf, 0.2, NA, 0.8, Inf, NA)))
    }, finally = endCluster())

多线程涉及更多的设置开销，因此只对非常大的栅格有意义and/or更复杂的计算（事实上，我很惊讶地注意到它并没有像在我下面测试的任何条件下的最佳选择——也许有更复杂的重新分类？）。

为了说明，我使用 OP 的设置绘制了微基准测试的结果，间隔高达 1000 万个单元格（每个单元格运行 10 次）如下：

最后一点，编译对任何测试尺寸都没有影响。

select 栅格数据有效范围的最快方法

Fastest way to select a valid range for raster data

r

r-raster