самый быстрый способ проведения арифметики c на больших растрах - PullRequest
Новая
       104

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

3 голосов
/ 10 июля 2020

У меня есть большое количество больших растров (глобальный размер, разрешение 250 м; около 1e10 ячеек с плавающей запятой) - имена файлов находятся в векторе deltaX.files. Я хочу добавить каждый из них к другому растру с именем X.tif. Поскольку эта операция может занять несколько дней, мне интересно, какой из них самый быстрый способ добавить растры, чтобы сделать это как можно быстрее.

Я могу придумать несколько методов, но я не уверен, какой из них наиболее эффективный вариант, или есть другой вариант лучше, чем любой из них.

Итак, у меня вопрос, есть ли способ оптимизировать или значительно ускорить арифметику c на больших растрах. Обратите внимание, что у меня есть графический процессор NVidia с поддержкой CUDA, поэтому очень приветствуются решения, которые могут распараллеливать это на графическом процессоре. Обратите внимание, что я использую Linux ystsem.

Некоторые примеры методов:

Обратите внимание на следующий блок кода, который нужно вставить перед каждым из них, чтобы определить сжатие выходного файла по умолчанию, выделение памяти и запуск параллельного кластера 

rasterOptions(chunksize = 1e10, maxmemory = 4e10)
f.opt = '-co COMPRESS=ZSTD -co PREDICTOR=2'
f.type = 'FLT4S'
beginCluster()

Опция (1) 

for (f in deltaX.files) {
  s = stack('X.tif', f)
  calc(s, sum, filename = paste0('new_', f), datatype = f.type, options = f.opt)
}

Вариант (2) 

X = raster('X.tif')
for (f in deltaX.files) {
  dX = raster(f)
  overlay(X, dX, fun=sum, filename = paste0('new_', f), datatype = f.type, options = f.opt)
}

Вариант (3) 

X = raster('X.tif')
for (f in deltaX.files) {
  dX = raster(f)
  Y = X + dX
  writeRaster(Y, filename = paste0('new_', f), datatype = f.type, options = f.opt)
}

Вариант (4): используйте gdal_cal c .py вместо R 

for (f in deltaX.files) {
  system(cmd)
  cmd = paste0("gdal_calc.py -A X.tif ", "-B ", f, " --outfile=", 'temp.tif', ' --calc="A+B"')
  system(cmd)
  system(paste('gdal_translate -ot Float32', f.opt, 'temp.tif', paste0('new_', f)))
  system('rm temp.tif')
}

Обратите внимание, что у меня были проблемы с тем, чтобы эта последняя версия успешно создавала полностью сжатые выходные файлы, поэтому также требуется дополнительный шаг использования gdal_translate для каждого файла для его сжатия. Однако на нескольких тестовых прогонах кажется, что значения искажены, поэтому меня больше всего интересует решение R, а не использование gdal_calc.py.

Некоторые фиктивные данные для воспроизводимости 

X = raster(vals = rnorm(65000 * 160000), ncol = 160000, nrow = 65000)
writeRaster(X, 'X.tif',  datatype = f.type, options = f.opt)
for (i in 1:10) {
  dX = raster(vals = rnorm(65000 * 160000), ncol = 160000, nrow = 65000)
  writeRaster(X, paste0('dX', i, '.tif'),  datatype = f.type, options = f.opt)
}
deltaX.files = paste0('dX', 1:10, '.tif')

1 Ответ

2 голосов
/ 10 июля 2020

Я бы посоветовал использовать terra (новый пакет, призванный заменить raster ---- он проще и быстрее). Теперь он доступен из CRAN, но для передового опыта вы можете установить с github

Вероятно, лучший подход - 

library(terra)
r <- rast(c('X.tif')
for (f in deltaX.files) {
    s <- rast(f)
    x <- c(r, s)
    y <- app(x, sum, filename=paste0('new_', f), datatype="INT2S",
           wopt=list(gdal="COMPRESS=LZW") )
 }

возможно ниже немного быстрее; но загвоздка в том, что у него нет аргумента имени файла. Но вы можете обойти это 

library(terra)
r <- rast(c('X.tif')
for (f in deltaX.files) {
    s <- rast(f)
    x <- r + s
    tempfile <- sources(x)$source[1]
    file.rename(tempfile, paste0('new_', f))
 }

В качестве альтернативы, за один шаг (это создаст один огромный файл - вероятно, нежелательно): 

r <- rast(c('X.tif')
s <- rast(deltaX.files)
# combine them as separate sub-datasets
x <- sds(r, s)
y <- sum(x, filename="file.tif")

Или вот так (быстро , но он попадает во временный файл, который вы можете переименовать, когда закончите, но вы не можете установить все параметры записи) 

 z <- r + s

Пока нет поддержки GPU ...

...