高效处理 r 中的大型空间数据答案

【问题标题】：Efficiently processing large spatial data in r高效处理 r 中的大型空间数据
【发布时间】：2021-01-12 23:57:43
【问题描述】：

我正在努力在以下代码中进行高效的数据处理。该代码产生了我想要的结果，一个包含两个变量的数据框，但它非常慢，因为空间数据集非常大，而且我是一个 R/编程初学者。该代码已经运行了大约 5 个小时，并且只循环了 488 个中的大约 80 个。

以下代码旨在确定在 488 个防火区中的每一个中被模拟野火横穿的建筑物的最大数量。火域是空间上连续的多边形，每个大约 100-300K 英亩。有超过一百万个模拟野火存储在 25 个不同的 shapefile 中，这些文件对应于称为火灾发生区 (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 GB

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

模拟火灾存储在与编号火灾发生区域 (FOA) 相关的 25 个独特形状文件中。

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 geospatial sf memory-efficient

【解决方案1】：

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 to subset a shapefile by a lat-long bounding box?。

降低%>% st_set_precision(0.05) %>% 的精度会减少运行时间吗？

祝你好运！

【讨论】：