эффективное использование растровых функций в r

Question

У меня есть 500+ очков в объекте SpatialPointsDataFrame; У меня есть 1,7 ГБ (200 000 строк х 200 000 колос) raster объект. Я хочу иметь табуляцию значений растровых ячеек в буфере вокруг каждой из 500 + точек.

Мне удалось добиться этого с помощью приведенного ниже кода (у меня было много вдохновения from here.). Тем не менее, он работает медленно, и я хотел бы, чтобы он работал быстрее. Фактически он работает нормально для буферов с «малой» шириной, скажем, 5 км до даже 15 км (~ 1 млн. Ячеек), но он становится очень медленным, когда буфер увеличивается до 100 км (~ 42 млн. Клеток).

Я мог бы легко улучшить петлю ниже, используя что-то из семейства apply и/или параллельного цикла. Но мое подозрение в том, что оно медленное, потому что пакет raster пишет 400 Мбайт + временные файлы для каждого взаимодействия цикла.

# packages 
library(rgeos) 
library(raster) 
library(rgdal) 

myPoints = readOGR(points_path, 'myLayer') 
myRaster = raster(raster_path) 

myFunction = function(polygon_obj, raster_obj) { 
    # this function return a tabulation of the values of raster cells 
    # inside a polygon (buffer) 

    # crop to extent of polygon 
    clip1 = crop(raster_obj, extent(polygon_obj)) 

    # crops to polygon edge & converts to raster 
    clip2 = rasterize(polygon_obj, clip1, mask = TRUE) 

    # much faster than extract 
    ext = getValues(clip2) 

    # tabulates the values of the raster in the polygon 
    tab = table(ext) 

    return(tab) 
} 

# loop over the points 
ids = unique(myPoints$ID) 
for (id in ids) { 

    # select point 
    myPoint = myPoints[myPoints$ID == id, ] 

    # create buffer 
    myPolygon = gBuffer(spgeom = myPoint, byid = FALSE, width = myWidth) 

    # extract the data I want (projections, etc are fine) 
    tab = myFunction(myPolygon, myRaster) 

    # do stuff with tab ... 
} 

Мои вопросы:

1. Правильно ли я частично виноваты операции письма? Если мне удастся избежать всех этих операций записи, будет ли этот код работать быстрее? У меня есть доступ к машине с 32 ГБ ОЗУ - так что я уверен, что я могу загрузить raster в память и не нужно писать временные файлы?

2. Что еще я мог сделать для повышения эффективности этого кода?

Answer 1

Я думаю, вы должны подойти к нему, как этот

library(raster) 
library(rgdal) 
myPoints <- readOGR(points_path, 'myLayer') 
myRaster <- raster(raster_path) 
e <- extract(myRaster, myPoints, buffer=myWidth) 

А потом что-то вроде

etab <- sapply(e, table) 

трудно ответить на ваш вопрос # 1, так как мы не достаточно знаем о вашем данных (мы не знаем, сколько ячеек покрыто буфером «100 км»). Но вы можете установить параметры, когда писать в файл с помощью функции rasterOptions. Вы заметили, что getValues быстрее, чем извлечение, на основе сообщения, на которое вы ссылаетесь, но я думаю, что это неправильно или, по крайней мере, не очень важно. Сочетание crop, rasterize и getValues должно иметь такую ​​же производительность, как и extract (что делает практически то, что под капотом). Если вы все равно идете по этому маршруту, вы должны передать пустой RasterLayer, созданный raster(myRaster) для более быстрого обрезки.