Эффективный способ извлечь область растра в R

Question

Я пытаюсь извлечь область ниже 5000 футов для области из 5 штатов, используя NED (National Elevation Dataset) с разрешением 30 метров (1 шт.), По штатам. У меня есть шейп-файл состояний и большой (около 4 ГБ) растр NED. Я использую пакет velox, который значительно ускоряет процесс, но я чувствую, что должен быть более быстрый способ сделать это (это все еще занимает много времени).

С пакетом velox я делаю растр в объект Velox, а затем извлечь клетки с помощью функции извлечения Velox. Поскольку я знаю разрешение каждой ячейки, я считаю количество ячеек ниже 5000 футов в каждом состоянии и умножаю на площадь каждой ячейки.

DEM_velox <- velox(NED)
DEM_cells <- DEM$extract(DEM_velox,sp=state_shapefile)

Нет ли лучшего и гораздо более быстрого способа сделай это? Я знаю, что есть функция area из растрового пакета и функция gArea из пакета rgeos, но, насколько я могу судить, они предназначены для полигонов (шейп-файлов), а не для растров. Любая помощь в этом была бы отличной.

Answer 1

Это пример решения вашей проблемы для Калифорнии с использованием al oop с набором более мелких растров (всего 69), которые охватывают весь штат и избегают extract.

растровых DEM. от шлюза пространственных данных NRCS.

Функция extract была исключена через mask каждого растра, который устанавливает высоты за пределами государственной границы равными NA, прежде чем подсчитать количество ячеек ниже критической отметки, используя cellStats вместе с функцией sum. Это суммирует количество ячеек ниже 5000 ', потому что растр преобразуется в значения TRUE / FALSE через ras_masked <= crit_elev до суммирования. Отдельный тест показал, что подход extract был в 5-6 раз медленнее, чем подход mask.

Время выполнения для l oop составляло 400 секунд.

library(raster)
crit_elev <- 5000 / 3.2808 #critical elevation in meters
demDir <- yourDEMdirectoryHERE
CADir <- yourStateBoundaryDirectoryHERE
list.files(demDir)
fnames <-list.files(demDir, full.names = TRUE)[grepl('*.tif$', list.files(demDir))]
length(fnames) #69 rasters that cover CA but some do spill over state boundary, so mask is necessary

#read in state boundary and project to DEM crs
list.files(CADir)
CA_TA.shp <- shapefile(file.path(CADir, 'california_CA_TA.shp'))
CA_UTM11N <- spTransform(CA_TA.shp, crs(raster(fnames[1])))

#loop to solve the problem for CA
elev_count <- vector(mode='integer', length = length(fnames)) #this is where we will keep track of how many cells in each dem are below the crit_elev
system.time(for (i in seq_along(fnames)) {
  ras_masked <- mask(raster(fnames[i]), CA_UTM11N)
  elev_count[i] <- cellStats(ras_masked <= crit_elev, sum)
  print(i)
})
#user  system elapsed 
#175.76  224.67  400.70 
#count cells below 5000 ft, multiply by resolution (30 m) and divide by 10,000 to get hectares
sum(elev_count) * 30 * 30 / 10000 
#[1] 19247149