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

Я борюсь с эффективной обработкой данных в следующем коде. Код дает желаемый результат - фрейм данных с двумя переменными, но он невероятно медленный, потому что наборы пространственных данных очень большие и я начинаю программировать на языке R /. Код работает около 5 часов и прошел всего около 80 из 488.

Следующий код предназначен для определения максимального количества зданий, пересекаемых смоделированными лесными пожарами в каждом из 488 очагов пожара. Пожарные навесы представляют собой пространственно смежные многоугольники, каждый площадью около 100–300 тыс. Акров. Существует более миллиона смоделированных лесных пожаров, хранящихся в 25 различных шейп-файлах, соответствующих географическим единицам планирования, называемым зонами возникновения пожаров (FOA), которые в среднем составляют 4-6 миллионов акров. Пожарные бассейны представляют собой небольшие полигоны на уровне сообществ, и их может быть несколько десятков в пределах более крупных региональных FOA.

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

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

  • Единицей анализа являются пожарные бассейны. Имеется 488 пожарных навесов (полигонов) *
PRJ <- st_crs(st_read("./Inputs/QWRAOriginal/FOA_perims/Original/foa401rfV3_Perims.shp"))

fshed <- st_read("./Inputs/Firesheds/QWRA_Firesheds.shp") %>% st_transform(crs=PRJ)

Это большой набор данных Microsoft Building Footprints. Это 1,8 ГБ

MBF_PNW <- st_read("./Inputs/MBF/MBF_PNW.shp")

Смоделированные пожары хранятся в 25 уникальных шейп-файлах, связанных с пронумерованными зонами возгорания (FOA).

FOA_table <- data.frame(FOA=c(401:412, 413, 413, 413,414:423),
                        File=list.files("./Inputs/QWRAOriginal/FOA_perims/Original/", pattern = ".shp$", full.names = TRUE))

#These polygons define the FOAs 
Bndry <- st_read("./Inputs/QWRAOriginal/FOA_boundaries.shp")

# empty output for storing loop products
outputs_table <- data.frame(Fireshed_N = NULL, Max_bldgs = NULL)

В следующем цикле: для каждого очага пожара определите любой FOA, который он пересекает, чтобы определить, какие моделируемые пожары следует считывать. В большинстве случаев это будет только один FOA, но в некоторых случаях очаги пожара могут быть расположены на границе между двумя или более FOA. Затем, определив соответствующие имитированные пожары, найдите все пожары, которые фактически пересекают очаг, и подсчитайте, сколько зданий было пересечено каждым пожаром. Меня интересует только максимальное количество зданий, пересекаемых одним пожаром, которое я затем сохраняю во фрейме выходных данных вместе с именем очага пожара.

for (i in 1:488){
  shp <- fshed[i,] #identify the polygon that represents fireshed 'i'

  foa_int <- Bndry %>% filter(st_intersects(., shp, sparse=FALSE)) %>% st_drop_geometry() %>% dplyr::select(FOA_number) %>% deframe() 
  
  tbl <- FOA_table %>% filter(FOA %in% c(unique(foa_int)))
  
  for (k in nrow(tbl)){
    fires <- st_read(tbl$File[k]) %>% filter(st_intersects(., shp, sparse=FALSE)) %>% 
      st_buffer(0) %>% st_set_precision(0.05) %>% st_make_valid() %>%  st_intersection(shp %>% st_make_valid) %>% 
      mutate(bldg_count=lengths(st_intersects(.,MBF_PNW))) %>% st_drop_geometry() 
    
    outputs_table <- rbind(outputs_table, data.frame(Fireshed_N=shp$Frshd_N, Max_bldgs=max(fires$bldg_count))) 
  }
  if(any(i==c(61,122,183,244,305,366,427,488))){
    write.csv(outputs_table, "./Outputs/Fshed_MBF.csv")* periodically save the output*
  }
  cat('Completed Fireshed',i) # add \n to the code
}

r sf geospatial memory-efficient
person Andy    schedule 12.01.2021    source источник

Ответы (1)


arrow_upward
0
arrow_downward

https://github.com/rstudio/profvis

Чтобы запустить код с профилированием, оберните выражение в profvis (). В результате появится интерактивный визуализатор профиля. Образец кода:

library(profvis)
library(ggplot2)

p <- profvis({
  g <- ggplot(diamonds, aes(carat, price)) + geom_point(size = 1, alpha = 0.2)
  print(g)
})


# View it with:
p
# or print(p)

Я бы изменил цикл for так, чтобы он запускался один раз, и профилировал функцию, чтобы выявить узкие места. (Полагаю, это будет эта команда st_intersects(.,MBF_PNW))

Можете ли вы предварительно вычислить какие-либо данные?

Можете ли вы фильтровать или подгруппировать свои данные, в частности "./Inputs/MBF/MBF_PNW.shp". См. R / GIS: How задать подмножество шейп-файла ограничивающей рамкой по широте и долготе?.

Приведет ли уменьшение точности %>% st_set_precision(0.05) %>% к сокращению времени выполнения?

Удачи!

person M.Viking    schedule 13.01.2021