大型栅格数据的矩阵索引匹配

Question

我有一个尺寸为 32251*51333 的大型栅格数据 (X)。 X 的值是另一个数组 (Y) 的重复，其大小为 3*10^6。 现在我想通过将 X 的值与 Y 的每个值进行匹配来更改 X 的值，例如我可以这样编程，

for (i in 1:length(Y)){
 X[X==Y[i]] = Z[i]   #Z is just another array with the same size as Y
}

问题在于，首先匹配X[X==Y[i]] = Z[i] 的索引不起作用，因为X 太大。几分钟后程序停止并给出错误"Error: cannot allocate vector of size 6.2 Gb". 其次，从 1 到长度（Y）的循环，即使 Y 的大小为 10^6，也可能需要“永远”才能完成。

我想到的一种方法是将 X 分成小块，然后为每个块进行索引匹配。但我觉得这仍然需要很多时间。

有没有更好的方法来实现上述目标？

第一次更新：

感谢@Lyngbakr 提供的示例，我将进一步阐述这个问题。因为我正在使用的栅格非常大（32251*51333），所以似乎无法上传。 @Lyngbakr 给出的示例与我想要的非常相似，只是创建的栅格太小。现在按照这个例子，我通过生成一个尺寸为 3000*2700 的更大的栅格来运行两个测试。请参阅下面的代码。

#Method 1: Use subs
start_time <- Sys.time()
Y <- 1:9
Z <- 91:99
X <- raster(matrix(rep(Y, 3), nrow=3000,ncol = 2700))
df <- data.frame(Y, Z)
X <- subs(X, df)
end_time <- Sys.time()
end_time - start_time
#Time difference of 2.248908 mins

#Method 2: Use for loop
start_time <- Sys.time()
Y <- 1:9
Z <- 91:99
X <- raster(matrix(rep(Y, 3), nrow=3000,ncol = 2700))
for (i in 1:length(Y)){
  X[X==Y[i]]=Z[i] #this indexing of R seems not efficient if X becomes large
}
end_time <- Sys.time()
end_time - start_time
#Time difference of 10.22717 secs

如您所见，一个简单的 for 循环甚至比 subs 函数更有效。请记住，示例中显示的栅格仍然小于我使用的栅格（大约小 100 个数量级）。此外，示例中的数组 Y 非常小。现在的问题可能是，如何加快方法 2，这只是一个简单的 for 循环？

Answer 1

您正在寻找subs 函数。我不知道它是否适用于大型栅格，但您可以尝试以下方法。

我加载了raster 包并创建了一些虚拟数据。 （如果您在问题中提供数据，这将是真的很有帮助。）然后，我绘制结果。

# Load library
library(raster)
#> Loading required package: sp

# Z holds values that will replace Y
Y <- 1:9
Z <- 91:99

# Create dummy raster
X <- raster(matrix(rep(Y, 3), ncol = 9))

# Examine raster
plot(X)

如您所见，X 只是一堆拼凑在一起的Y 向量。接下来，我将Y 和Z 绑定到一个数据框df。

# Combine y & z into a data frame
df <- data.frame(Y, Z)

最后，我使用subs 将Y 值替换为Z 值。

# Substitute Z for Y in X
X <- subs(X, df)

快速查看栅格显示值已被正确替换。

# Examine raster
plot(X)

^{由reprex package (v0.2.1.9000) 于 2019 年 6 月 25 日创建}

---

更新

Rcpp在性能成为问题时非常有用。下面，我比较了三种方法：

1. 在 R 中循环（来自问题）
2. 使用光栅包中的subs
3. 使用 Rcpp 在 C++ 中循环

顺便说一句，Sys.time() 不是检查性能的好方法，所以我建议改为使用microbenchmark。

# Load library
library(raster)

# Define vectors and raster
Y <- 1:9
Z <- 91:99
X <- raster(matrix(rep(Y, 3), nrow = 3000, ncol = 2700))

method_1 是subs 函数。

# Using subs function
method_1 <- function(){
  df <- data.frame(Y, Z)
  X <- subs(X, df)
}

method_2 是您最初的循环方法。

# Using R loop
method_2 <- function(){
  for (i in 1:length(Y)){
    X[X==Y[i]]=Z[i] 
  }
  X
}

method_3 是用 C++ 实现的循环方法。

# Using Rcpp loops
src <-
"Rcpp::NumericMatrix subs_cpp(Rcpp::NumericMatrix X, Rcpp::NumericVector Y, Rcpp::NumericVector Z){
  for(int i = 0; i < Y.length(); ++i){
    for(int j = 0; j < X.ncol(); ++j){
      for(int k = 0; k < X.nrow(); ++k){
        if(X(k, j) == Y(i)){
          X(k, j) = Z(i);
        }
      }
    }
  }  

  return X;
}"

Rcpp::cppFunction(src)

method_3 <- function(){
  subs_cpp(as.matrix(X), Y, Z)
}

我在这里对方法进行基准测试。

# Run benchmarking
microbenchmark::microbenchmark(method_1(), method_2(), method_3(), times = 10)

# Unit: milliseconds
#       expr        min         lq       mean     median         uq       max neval
# method_1() 16861.5447 17737.2124 19321.5674 18628.8573 20117.0159 25506.208    10
# method_2()   671.2223   677.6029  1111.3935   738.6216  1657.0542  2163.137    10
# method_3()   316.9810   319.1484   481.3548   320.2337   326.7133  1477.454    10

如您所见，Rcpp 方法是迄今为止最快的。

您还可以比较输出以确保它们使用较小的栅格产生相同的结果。

# Examine all three outputs with smaller raster
X <- raster(matrix(rep(Y, 3), ncol = 9))

plot(method_1(), main = "Method 1")
plot(method_2(), main = "Method 2")
plot(raster(method_3()), main = "Method 3") # Needs to converted into a raster

而且它们看起来都很相似。请注意，对于第三种方法，需要将结果从矩阵转换回栅格。