You can use spatstat
to do this quite efficiently if it is OK to treat your
lat,lon data as flat coordinates. So we assume you don’t cross the date line
(+/-180 degrees longitude) and that 1 degree latitude is equal to 1 degree
longitude. That is only correct around equator. Away from the equator you
should really multiply one of the coordinates with an appropriate factor or
you should ideally project your data to a flat coordinate system.
With randomly generated data of 16 stations and almost 3 million depth
measurements the code finding nearest neighbours below takes less than half
a second on my laptop.
library(spatstat)
stations <- data.frame(stationID = 1:16,
Latitude = runif(16, 0, 10),
Longitude = runif(16, 0, 10))
Wstations <- bounding.box.xy(stations$Longitude, stations$Latitude)
Xstations <- ppp(stations$Longitude, stations$Latitude, window = Wstations)
N <- 2943485
depths <- data.frame(Longitude = runif(N, 0, 10),
Latitude = runif(N, 0, 10),
depth = runif(N, 0, 1000))
Wdepths <- bounding.box.xy(depths$Longitude, depths$Latitude)
Xdepths <- ppp(depths$Longitude, depths$Latitude, window = Wdepths)
id <- nncross(Xstations, Xdepths, what = "which")
stations$depths <- depths$depth[id]
head(stations)
#> stationID Latitude Longitude depths
#> 1 1 5.7992147 5.6716015 435.1266
#> 2 2 9.2218643 6.0519959 154.5833
#> 3 3 7.6444158 3.2228619 963.7626
#> 4 4 3.8993755 0.9428149 204.9189
#> 5 5 7.9673171 0.9801396 933.5696
#> 6 6 0.9453616 8.0829834 603.0325