（空间）找到点X米内所有点的有效方法？
Posted
技术标签:
【中文标题】（空间）找到点X米内所有点的有效方法？【英文标题】：(Spatial) Efficient way of finding all points within X meters of a point? 【发布时间】：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
【问题讨论】：
你能举一个最小的例子让我们试试吗？见***.com/questions/5963269/… 抱歉，我认为我提供的示例代码应该足够了？我发布的示例不符合可复制示例的标准怎么办？ @Tim_K 最后我很好奇，我实现了一个集成的 sf + data.table 可能的解决方案。您可能对下面的更新答案感兴趣。 你应该考虑看看这个帖子：gis.stackexchange.com/questions/255671/…;我遇到了同样的问题，并用近似值和data.table 子集解决了它，它也可以很容易地并行运行。我不确定这是否是最快的方法，但对于 9*10^6，单核大约需要 80 小时，2 核需要 40 小时等等。
nilsole 该帖子有助于思考问题。提出的解决方案是在进行多边形点计算之前使用平方子集进行预过滤。类似于下面@lbusett 的答案，但是，子集是在每个单独的点上完成的，而不是将整个平面雕刻成一个 nxn 网格
【参考方案1】：
我有两种选择，一种看起来更快，另一种则不然。不幸的是，更快的方法可能不适用于并行化，因此它可能无济于事。
library(sf)
nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 32618)
# create points
pts <- st_centroid(nc)
dis <- 50000
result <- list()
你的方法
system.time(
for (i in 1:nrow(pts))
b <- st_buffer(pts[i,], dist = dis)
result[[i]] <- st_intersects(b, nc)[[1]]
)
较慢的选择
system.time(
for (i in 1:nrow(pts))
b <- as.vector(st_distance(pts[i,], pts))
result[[i]] <- which(b <= dis)
)
对于较小的数据集，没有循环：
x <- st_distance(pts)
res <- apply(x, 1, function(i) which(i < dis))
更快的替代方案（不明显如何并行执行），并且可能是不公平的比较，因为我们自己不进行循环
library(spdep)
pts2 <- st_coordinates(pts)
system.time(x <- dnearneigh(pts2, 0, dis))
我会首先得到一个包含指示邻居的索引的列表，然后提取属性（应该很快）
【讨论】：
根据您的回答，我能够找到这篇博文，进一步讨论了同一主题：cran.r-project.org/web/packages/spdep/vignettes/nb_sf.html 与以上可以在停留在 sf 内时应用，例如 x 【参考方案2】：根据 RobertH 的回答，在这个特定示例中使用 sf::st_coordinates 提取坐标要快一些。
library(sf)
library(spdep)
nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 32618)
# create points
pts <- st_centroid(nc)
dis <- 50000
# quickest solution:
x <- spdep::dnearneigh(sf::st_coordinates(pts), 0, dis)
微基准测试：
my_method <- function(pts)
result <- list()
for (i in 1:nrow(pts))
b <- st_buffer(pts[i,], dist = dis)
result[[i]] <- st_intersects(b, nc)[[1]]
result
library(microbenchmark)
microbenchmark(
my_method(pts),
dnearneigh(as(pts, 'Spatial'), 0, dis),
dnearneigh(st_coordinates(pts), 0, dis)
)
Unit: microseconds
expr min lq mean median uq max neval
my_method(pts) 422807.146 427434.3450 435974.4320 429862.8705 434968.3975 596832.271 100
dnearneigh(as(pts, "Spatial"), 0, dis) 3727.221 3939.8540 4155.3094 4112.8200 4221.9525 7592.739 100
dnearneigh(st_coordinates(pts), 0, dis) 394.323 409.5275 447.1614 430.4285 484.0335 611.970 100
检查等价性：
x <- dnearneigh(as(pts, 'Spatial'), 0, dis)
y <- dnearneigh(st_coordinates(pts), 0, dis)
all.equal(x,y, check.attributes = F)
[1] TRUE
【讨论】：
as(pts, 'Spatial') 将 sf 对象转换为 sp. 中定义的 Spatial* 对象。它不是spdep 的一部分。 dnearneigh 接受坐标矩阵的空间对象。提取坐标更快，但两种方法都很快，您只需为整个数据集执行一次，因此差异不应该那么重要。 （它应该或多或少地线性缩放 --- 而距离计算没有）
你说得对。我调整了答案中的语言来解决这个问题。我上面的例子是针对这个用例的，不一定适用于一般情况。【参考方案3】：
使用sf 对象时，显式循环要执行的功能
相交之类的二元运算通常会适得其反（另请参阅
How can I speed up spatial operations in `dplyr::mutate()`?)
类似于您的方法（即缓冲和相交），但没有
显式的for 循环效果更好。
让我们看看它在一个相当大的 50000 点数据集上的表现：
library(sf)
library(spdep)
library(sf)
pts <- data.frame(x = runif(50000, 0, 100000),
y = runif(50000, 0, 100000))
pts <- sf::st_as_sf(pts, coords = c("x", "y"), remove = F)
pts_buf <- sf::st_buffer(pts, 5000)
coords <- sf::st_coordinates(pts)
microbenchmark::microbenchmark(
sf_int = int <- sf::st_intersects(pts_buf, pts),
spdep = x <- spdep::dnearneigh(coords, 0, 5000)
, times = 1)
#> Unit: seconds
#> expr min lq mean median uq max neval
#> sf_int 21.56186 21.56186 21.56186 21.56186 21.56186 21.56186 1
#> spdep 108.89683 108.89683 108.89683 108.89683 108.89683 108.89683 1
您可以在此处看到st_intersects 方法比
dnearneigh 一个。
很遗憾，这不太可能解决您的问题。看执行 我们得到的不同大小的数据集的时间：
subs <- c(1000, 3000, 5000, 10000, 15000, 30000, 50000)
times <- NULL
for (sub in subs[1:7])
pts_sub <- pts[1:sub,]
buf_sub <- pts_buf[1:sub,]
t0 <- Sys.time()
int <- sf::st_intersects(buf_sub, pts_sub)
times <- cbind(times, as.numeric(difftime(Sys.time() , t0, units = "secs")))
plot(subs, times)
times <- as.numeric(times)
reg <- lm(times~subs+I(subs^2))
summary(reg)
#>
#> Call:
#> lm(formula = times ~ subs + I(subs^2))
#>
#> Residuals:
#> 1 2 3 4 5 6 7
#> -0.16680 -0.02686 0.03808 0.21431 0.10824 -0.23193 0.06496
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.429e-01 1.371e-01 1.772 0.151
#> subs -2.388e-05 1.717e-05 -1.391 0.237
#> I(subs^2) 8.986e-09 3.317e-10 27.087 1.1e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.1908 on 4 degrees of freedom
#> Multiple R-squared: 0.9996, Adjusted R-squared: 0.9994
#> F-statistic: 5110 on 2 and 4 DF, p-value: 1.531e-07
在这里，我们看到时间和时间之间几乎完美的二次关系 点数（正如预期的那样）。在 10M 点子集上，假设 行为没有改变，你会得到：
predict(reg, newdata = data.frame(subs = 10E6))
#> 1
#> 898355.4
，对应10天左右，假设趋势不变
当进一步增加点数时（但同样会发生
dnearneigh...)
我的建议是将你的观点“拆分”成块，然后在一个 每个分割的基础。
例如，您可以在开始时对您的积分进行排序 x 轴，然后使用 data.table 轻松快速地提取缓冲区和点的子集与它们进行比较。
显然，“点”缓冲区需要大于“缓冲区”，根据
到比较距离。因此，例如，如果您使用 pts_buf 制作一个子集
[50000 - 55000] 中的质心，pts 的相应子集应包括
[49500 - 55500] 范围内的点。
通过将不同的子集分配给
foreach 或类似构造中的不同内核。
我什至不知道在这里使用空间对象/操作是否有益，因为一旦我们有了坐标，我们所需要的就是计算和子集欧几里德距离：我怀疑基于 data.table 的仔细编码的蛮力方法可能是也是一个可行的解决方案。
HTH！
更新
最后，我决定试一试，看看我们可以从这种方法中获得多少速度。这是一个可能的实现：
points_in_distance_parallel <- function(in_pts,
maxdist,
ncuts = 10)
require(doParallel)
require(foreach)
require(data.table)
require(sf)
# convert points to data.table and create a unique identifier
pts <- data.table(in_pts)
pts <- pts[, or_id := 1:dim(in_pts)[1]]
# divide the extent in quadrants in ncuts*ncuts quadrants and assign each
# point to a quadrant, then create the index over "xcut"
range_x <- range(pts$x)
limits_x <-(range_x[1] + (0:ncuts)*(range_x[2] - range_x[1])/ncuts)
range_y <- range(pts$y)
limits_y <- range_y[1] + (0:ncuts)*(range_y[2] - range_y[1])/ncuts
pts[, `:=`(xcut = as.integer(cut(x, ncuts, labels = 1:ncuts)),
ycut = as.integer(cut(y, ncuts, labels = 1:ncuts)))] %>%
setkey(xcut, ycut)
results <- list()
cl <- parallel::makeCluster(parallel::detectCores() - 2, type =
ifelse(.Platform$OS.type != "windows", "FORK",
"PSOCK"))
doParallel::registerDoParallel(cl)
# start cycling over quadrants
out <- foreach(cutx = seq_len(ncuts)), .packages = c("sf", "data.table")) %dopar%
count <- 0
# get the points included in a x-slice extended by `dist`, and build
# an index over y
min_x_comp <- ifelse(cutx == 1, limits_x[cutx], (limits_x[cutx] - maxdist))
max_x_comp <- ifelse(cutx == ncuts,
limits_x[cutx + 1],
(limits_x[cutx + 1] + maxdist))
subpts_x <- pts[x >= min_x_comp & x < max_x_comp] %>%
setkey(y)
for (cuty in seq_len(pts$ycut))
count <- count + 1
# subset over subpts_x to find the final set of points needed for the
# comparisons
min_y_comp <- ifelse(cuty == 1,
limits_y[cuty],
(limits_y[cuty] - maxdist))
max_y_comp <- ifelse(cuty == ncuts,
limits_y[cuty + 1],
(limits_y[cuty + 1] + maxdist))
subpts_comp <- subpts_x[y >= min_y_comp & y < max_y_comp]
# subset over subpts_comp to get the points included in a x/y chunk,
# which "neighbours" we want to find. Then buffer them.
subpts_buf <- subpts_comp[ycut == cuty & xcut == cutx] %>%
sf::st_as_sf() %>%
st_buffer(maxdist)
# retransform to sf since data.tables lost the geometric attrributes
subpts_comp <- sf::st_as_sf(subpts_comp)
# compute the intersection and save results in a element of "results".
# For each point, save its "or_id" and the "or_ids" of the points within "dist"
inters <- sf::st_intersects(subpts_buf, subpts_comp)
# save results
results[[count]] <- data.table(
id = subpts_buf$or_id,
int_ids = lapply(inters, FUN = function(x) subpts_comp$or_id[x]))
return(data.table::rbindlist(results))
parallel::stopCluster(cl)
data.table::rbindlist(out)
函数将点sf对象、目标距离和数字作为输入
用于划分象限范围的“切割”，并在输出中提供
一个数据框，对于每个原始点，其中的点的“ID”
maxdist 在int_ids 列表列中报告。
在具有不同数量的均匀分布点的测试数据集上，
和maxdist 的两个值我得到了这样的结果（“并行”运行是使用 6 个内核完成的）：
因此，我们在“串行”实现上已经速度提高了 5-6 倍，并且由于超过 6 个内核的并行化，又提高了 5 倍。 虽然这里显示的时间只是指示性的，并且与 我们构建的特定测试数据集（在分布不太均匀的数据集上，我预计速度会降低）我认为这非常好。
HTH！
PS：更全面的分析可以在这里找到：
https://lbusettspatialr.blogspot.it/2018/02/speeding-up-spatial-analyses-by.html
【讨论】：
出于文档目的，我认为您答案顶部的 SO 问题中的这条评论似乎相关：“如果步骤涉及二进制逻辑谓词（如 st_intersects、st_crosses 等，请避免按行操作） .) 因为你失去了空间索引效率提升"以上是关于（空间）找到点X米内所有点的有效方法？的主要内容，如果未能解决你的问题，请参考以下文章