При работе с объектами sf явный цикл по функциям для выполнения бинарных операций, таких как пересечения, обычно контрпродуктивен (см. Также Как я могу ускорить пространственные операции в` 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 в 5 раз быстрее, чем метод 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
Здесь мы видим почти идеальную квадратичную зависимость между временем и количеством точек (как и следовало ожидать). На подмножестве из 10 миллионов точек, предполагая, что поведение не изменится, вы получите:
predict(reg, newdata = data.frame(subs = 10E6))
#> 1
#> 898355.4
, что соответствует примерно 10 дням, если предположить, что тренд будет постоянным при дальнейшем увеличении количества точек (но то же самое произойдет и для _9 _...)
Я предлагаю «разбить» ваши баллы на куски, а затем работать над ними.
Например, вы можете расположить свои точки в начале по оси 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, целевое расстояние и количество «разрезов», которые используются для разделения экстента на квадранты. и предоставляет на выходе кадр данных, в котором для каждой исходной точки "идентификаторы" точек в пределах maxdist сообщаются в столбце списка int_ids.
На тестовом наборе данных с различным количеством равномерно распределенных точек и двумя значениями maxdist я получил такие результаты («параллельный» прогон выполняется с использованием 6 ядер):
Итак, здесь мы получаем 5-6-кратное увеличение скорости уже при «последовательной» реализации и еще 5-кратное благодаря распараллеливанию по 6 ядрам. Хотя время, показанное здесь, является лишь ориентировочным и связано с конкретным набором тестовых данных, который мы создали (на менее равномерно распределенном наборе данных я бы ожидал более низкого улучшения скорости), я думаю, что это неплохо.
HTH!
PS: более подробный анализ можно найти здесь:
https://lbusettspatialr.blogspot.it/2018/02/speeding-up-spatial-analyses-by.html
person
lbusett
schedule
09.02.2018
data.tableподмножества, которые также можно легко запустить параллельно. Я не уверен, что это самый быстрый способ сделать это, но для 9 * 10 ^ 6 требуется около 80 часов на одно ядро, 40 часов на 2 ядра и так далее. - person nilsole schedule 22.02.2018