【发布时间】：2018-07-16 22:40:05
【问题描述】：
我有一个大型空间数据集（12M 行）。几何图形是地图上的点。对于数据集中的每一行，我想找到该点 500 米范围内的所有点。
在 r 中，使用 sf，我一直试图通过并行循环每一行并运行 st_buffer 和 st_intersects 来做到这一点，然后将结果以键值格式保存为列表（键是原点，值是邻居）。
问题是数据集太大。即使并行处理超过 60 个内核，操作也需要很长时间（>1 周并且通常会崩溃）。
这种蛮力方法的替代方法是什么？是否可以使用 sf 构建索引？也许将操作推送到外部数据库？
代表：
library(sf)
library(tidyverse)
library(parallel)
library(foreach)
# example data, convert to decimal:
nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% st_transform(32618)
# expand the data a a bit to make the example more interesting:
nc <- rbind(nc,nc,nc)
nc <- nc %>% mutate(Id = row_number())
## can run in parallel if desired:
# num_cores <- parallel::detectCores()-2
# cl <- makeSOCKcluster(num_cores)
# registerDoSNOW(cl)
# or just run in sequence:
registerDoSEQ()
neighbors <- foreach(ii = 1:nrow(nc)
, .verbose = FALSE
, .errorhandling = "pass") %dopar% {
l = 500 # 500 meters
# isolate the row as the origin point:
row_interest <- filter(nc, row_number()==ii)
# create the buffer:
buffer <- row_interest %>% st_buffer(dist = l)
# extract the row numbers of the neighbors
comps_idx <- suppressMessages(st_intersects(buffer, nc))[[1]]
# get all the neighbors:
comps <- nc %>% filter(row_number() %in% comps_idx)
# remove the geometry:
comps <- comps %>% st_set_geometry(NULL)
# flow control in case there are no neibors:
if(nrow(comps)>0) {
comps$Origin_Key <- row_interest$Id
} else {
comps <- data_frame("lat" = NA_integer_,"lon" = NA_integer_, "bbl" = row_interest$bbl)
comps$Origin_Key <- row_interest$Id
}
return(comps)
}
closeAllConnections()
length(neighbors)==nrow(nc)
[1] TRUE
【问题讨论】：
-
你能举一个最小的例子让我们试试吗？见stackoverflow.com/questions/5963269/…
-
抱歉，我认为我提供的示例代码应该足够了？我发布的示例不符合可复制示例的标准怎么办？
-
@Tim_K 最后我很好奇，我实现了一个集成的 sf + data.table 可能的解决方案。您可能对下面的更新答案感兴趣。
-
你应该考虑看看这个帖子：gis.stackexchange.com/questions/255671/…;我遇到了同样的问题，并用近似值和
data.table子集解决了它，它也可以很容易地并行运行。我不确定这是否是最快的方法，但对于 9*10^6，单核大约需要 80 小时，2 核需要 40 小时等等。 -
nilsole 该帖子有助于思考问题。提出的解决方案是在进行多边形点计算之前使用平方子集进行预过滤。类似于下面@lbusett 的答案，但是，子集是在每个单独的点上完成的，而不是将整个平面雕刻成一个 nxn 网格