# Nearest point to 30 million points using R

I have two point datasets. One with 30 million records and the other with 1.5 million. For each of the 30 million points, I would like to calculate the nearest point(s) from the 1.5 million dataset to them.

I have the following, very simple, code that would do the job if I had a spare 6 months. Is there a quick and efficient alternative way to do this?

``````library(data.table)

DT1 <- data.table(DT1_id = 40000001:70000000, DT1_x =sample(450000:550000, size = 30000000,replace = TRUE),DT1_y =sample(250000:350000, size = 30000000,replace = TRUE) )
DT2 <- data.table(DT2_id = 8000001:9500000, DT2_x =sample(450000:550000, size = 1500000,replace = TRUE),DT2_y =sample(250000:350000, size = 1500000,replace = TRUE) )

DT4 <- as.data.table(NULL)

for (n in 1:nrow(DT1)){
DT3 <- cbind(DT1[rep(n,nrow(DT2))],DT2)
DT3[, DIST := sqrt((DT1_x - DT2_x)^2 + (DT1_y - DT2_y)^2)]
DT4 <- rbind(DT4,DT3[ , .SD[DIST==min(DIST)], by = DT1_id])}
``````

EDIT

Revised DT2 example dataset that more closely resembles my real data:

``````DT2 <- data.table(DT2_x =sample(450000:550000, size = 1200000,replace = TRUE),DT2_y =sample(250000:350000, size = 1200000,replace = TRUE))
DT2 <- rbind(DT2,DT2[sample(.N, 400000)])
DT2 <- data.table(cbind(DT2_id = sprintf("%s%0*d", "UID_", 7, 1:nrow(DT2)),DT2))
``````

Using a mix of nearest neighbour search (RANN) and parallelization (future/furrr)

``````library(RANN)
lbrary(future)
library(furrr)
plan(multisession) # launch local cluster

# if you get message about limits in globals size then increase it
options(future.globals.maxSize=600*1024*1024)

# split DT1 in chunks for parallelization
chunk_size <- 1000000
DT1[, chunk:=floor(.I/chunk_size)]
chunked_DT1 <- split(DT1, by="chunk", keep.by = FALSE)

# apply compute function on each chunk using parallelization
res <- furrr::future_map(chunked_DT1, function(chunk) {
library(data.table)
neighbours <- RANN::nn2(DT2[,-1], chunk[,-1], 1) # 1 stands for 1 nearest neighbours
chunk[, DT2_id:=DT2[as.vector(neighbours\$nn.idx)]\$DT2_id] # get DT2 IDs from chunck IDs
})

# bind all results in an augmented DT1
res <- rbindlist(res)
``````
• This is a nice solution. Besides the use of RANN, I like the parallelization addition. Commented Mar 31, 2021 at 23:40
• @Billy34 that's great - thank you! It works great with my real data, except in cases where there are multiple DT2 points that are closest to a DT1 point. Can the solution be amended to deal with this (not sure where to start doing it myself as these are all new packages for me!)? Commented Apr 1, 2021 at 10:11
• Do you mean two or more DT2 points that are at the same exact distance of the DT1 one ? If so you can increase the max number of search neighbours (here as you wanted the closest, I set it to 1). It means that we need to provide a sensible value to it. nn2 also returns the distance (see documentation). So we need to detect such cases ans return corresponding IDs. By the way if my answer is of interest to you, would you mind accept it ? (grey checkmark close to my answer) Commented Apr 1, 2021 at 16:02
• I have added a revised DT2 dataset that more closely resembles my real data. Basically some points in DT2 will have the same XYs. In those cases, if they are the nearest to a DT1 point, I would like however many of those points that are the closest to be returned (i.e. 1 to many). In other cases where there is a single DT2 point then it will be a 1 to 1 return. So I imagine the final res object will have more rows than the original DT1. Sorry I didn't include this in my original example. Also - increasing the number of neighbours to be returned to 2 causes an error Commented Apr 1, 2021 at 18:47
• @Chris Please open a new question rather than tagging on another question to this one.
– Aaron
Commented Apr 1, 2021 at 19:14

If your DT2 dataset has multiple points with the same XY coordinates and you wish all of these to be included in your output (if they are closest to a particular point in DT1) then running this after @Billy34's answer should produce what you're looking for:

``````library(data.table)

res <- res[DT2, on = c(DT2_id = "DT2_id"), nomatch = 0]

res <- DT2[res, on = c(DT2_x="DT2_x", DT2_y="DT2_y"), nomatch = 0, allow.cartesian=TRUE
][,.(DT1_id, DT1_x, DT1_y, DT2_id, DT2_x, DT2_y)]

setorder(res, DT1_id, DT2_id)
``````