GIS-LAB

Географические информационные системы и дистанционное зондирование

Сотни тысяч пикселей и скорость операций в R

Максим Дубинин, 05.12.2007

R и rgdal классные штуки, но с обработкой большого количества данных, а растрами это обычно так и есть, надо быть осторожным. Допустим есть 2 растра, второй получен на основе первого и от некоторых его значений надо избавиться (несчастная черная рамка вокруг значимой части первого растра). Само собой линейные размеры первого второго одинаковые.

Возьмем простую операцию – условие, например такое: если для определенного пикселя сумма значений во всех каналах первого растра = 32000, присвоить значение 0 соответствующему пикселю второго растра.

В R это будет выглядеть вот так:

for (i in 1:length(raster1)) {
if (sum(raster1[ind,]) == 32000) {
raster2[ind] = 0
}
}

Но работает дооолго, нереально долго, начинаем разбираться что именно работает долго и методом проб и ошибок выясняем, что это сумма. Сделаем тест попроще, изменив условие: даже если только в первом канале 0, то и такое условие сгодится, перепишем:

for (i in 1:length(raster1)) {
if (sum(raster1[ind,1]) == 0) {
raster2[ind] = 0
}
}

Очень быстро.

Мораль (одна из…) – надо мерять время выполнения операций и оптимизировать узкие места. Как? System.time
Попробуем, единственный недостаток, цикл придется “разложить” в строку. Сравним время выполнения операций (количество пикселей 266612):

>system.time(for (ind in 1:length(raster1)) {if (sum(raster1[ind,]) == 32000) {raster2[ind] = 0}})
elapsed =  715.47
>system.time(for (ind in 1:length(raster1)) {if (sum(raster1[ind,1]) == 0) { raster2[ind] = 0}})
elapsed = 29.01

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

Arcinfo Workstation справляется с подобной задачкой шутя за 5 секунд + кое какие расходы на конвертацию данных, однако до изящности и лаконичности ей далековато и плодить гриды ох как не хочется:

raster3 = con(sum(raster1c1,raster1c2,raster1c3,raster1c4,raster1c5,
raster1c6,raster1c7,raster1c8,raster1c9,raster1c10,
raster1c11,raster1c12,raster1c13,raster1c14,raster1c15) == 32000,0, raster2)

Комментарии (1) к статье “Сотни тысяч пикселей и скорость операций в R”

  1. sim says:

    Ну а совсем это убыстрить можно так, перед условием избавиться от остальных колонок в проверяемом массиве:

    subraster1 = raster1[,1]
    system.time(for (ind in 1:length(raster1)) {if (img2[ind] == 32000) {raster2[ind] = 0}})[3]

    elapsed
    1.95

Оставьте комментарий


(Геокруг)

Если Вы обнаружили на сайте ошибку, выберите фрагмент текста и нажмите Ctrl+Enter