Nearest point to 30 million points using R

Question

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))

Answer 1

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)

Answer 2

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)