сопоставление индекса матрицы для больших растровых данных

Question

У меня есть большие растровые данные (X) с размером 32251 * 51333. Значения X являются повторениями другого массива (Y), который имеет размер 3 * 10 ^ 6. Теперь я хочу изменить значения X, сопоставляя его с каждым значением Y, например, я могу запрограммировать так:

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 на маленькие порции, а затем выполнить индексное сопоставление для каждого порции. Но я чувствую, что это все равно займет много времени.

Есть ли лучший способ достичь вышеуказанной цели?

1-е обновление:

Благодаря примеру, предоставленному @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

Как видите, простой цикл для даже более эффективен, чем функция sub. Помните, что растр, показанный в примере, все еще меньше, чем тот, с которым я работаю (примерно на 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)

^{Создано в 2019-06-25 пакетом Представить (v0.2.1.9000)}

---

Обновление

Rcpp действительно полезно, когда производительность является проблемой.Ниже я сравниваю три метода:

1. Цикл в R (из вопроса)
2. Использование subs из растрового пакета
3. Цикл в C ++ с использованием Rcpp

Кстати, 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

И все они выглядят одинаково.Обратите внимание, что для третьего метода результат необходимо преобразовать обратно в растр из матрицы.