function for finding all points within X meters of a point
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) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This comment has been minimized.
FYI - really love this modification of the points in distance function and was trying to use it for an experiment for a class. I ran into some issues running it initial (extra parentheses in the foreach()) and a typo in the inner for loop. It seemed to work with these changes: