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)
Ну а совсем это убыстрить можно так, перед условием избавиться от остальных колонок в проверяемом массиве:
subraster1 = raster1[,1]
system.time(for (ind in 1:length(raster1)) {if (img2[ind] == 32000) {raster2[ind] = 0}})[3]
elapsed
1.95