I've got a dataframe with latitude and longitude, which looks like this:
x y set
61 -112
63 -113
61 -113
62 -111 point
61 -111
64 -120
I want to find the three closest points to the point that is marked as point in column set. Then, for these three closest points, I want to amend the column set to say closest. Like this:
x y set
61 -112 closest
65 -113
62 -113 closest
62 -111 point
62 -111 closest
64 -120
How can I do this?
dists <- geosphere::distHaversine(dat[dat$set=="point",c("y","x")], dat[,c("y","x")])
dists
# [1] 123339.4 151513.9 153862.4 0.0 111319.5 505814.4
dat$set[dat$set != "point" & rank(dists) < 5] <- "closest"
dat
# x y set
# 1 61 -112 closest
# 2 63 -113 closest
# 3 61 -113
# 4 62 -111 point
# 5 61 -111 closest
# 6 64 -120
The reason we use < 5 is that the own-distance (point to point) will be the closest (0), so we need ranks 2-4. This assumes there is one "point"; if there are more, you'll likely want outer (to produce a matrix of distances) and look at each row before populating $set.
I'm inferring latitude and longitude from the sp tag, so chose the Haversine distance calculation since it's fast, and the appearance of coarse coordinates does not suggest the requirement for sub-millimeter accuracy (i.e., Vincenty Ellipsoid formula). There are other distance calculations if needed.
Here is first another approach with geosphere (make a distance matrix with distm) and then I show how you can use the terra::nearby method (which works for both long/lat and planar coordinates).
m <- matrix(c(61, -112, 63, -113, 61, -113, 62, -111, 61, -111, 64, -120), ncol=2, byrow=TRUE)
# note that the order should be long/lat !!!
m <- m[, 2:1]
d <- geosphere::distm(m)
diag(d) <- NA
i <- order(d[4,])[1:3]
i
#[1] 5 1 2
m[i,]
# [,1] [,2]
#[1,] -111 61
#[2,] -112 61
#[3,] -113 63
Now with terra. The below gets the nearest 3 neighbors for all points.
library(terra)
v <- vect(m, crs="+proj=lonlat")
nearby(v, k=3)
# id k1 k2 k3
#1 1 3 5 4
#2 2 4 3 1
#3 3 1 5 4
#4 4 5 1 2
#5 5 1 3 4
#6 6 2 3 4
With terra version 1.3.15 (currently the development version) you can also do
nearby(v[4,], v, k=4)
# id k1 k2 k3 k4
#[1,] 1 4 5 1 2
Taking k=4 neighbors as the first one is the point itself.
To get the development version, do
install.packages('terra', repos='https://rspatial.r-universe.dev')