Find closest points (lat / lon) from one data set to a second data set

Question

I have two data sets, A and B, which give locations of different points in the UK as such:

A = data.frame(reference = c(C, D, E), latitude = c(55.32043, 55.59062, 55.60859), longitude = c(-2.3954998, -2.0650243, -2.0650542))

B = data.frame(reference = c(C, D, E), latitude = c(55.15858, 55.60859, 55.59062), longitude = c(-2.4252843, -2.0650542, -2.0650243))

A has 400 rows and B has 1800 rows.
For all the rows in A, I would like to find the shortest distance in kilometers between a point in A and each of the three closest points in B, as well as the reference and coordinates in lat and long of these points in B.

I tried using this post

R - Finding closest neighboring point and number of neighbors within a given radius, coordinates lat-long

However, even when I follow all the instructions, mainly using the command distm from the package geosphere, the distance comes up in a unit that can't possibly be kilometers. I don't see what to change in the code, especially since I am not familiar at all with the geo packages.

Apps

Answer 1

Here is solution using a single loop and vectorizing the distance calculation (converted to km).
The code is using base R's rank function to order/sort the list of calculated distances.
The indexes and the calculated distances of the 3 shortest values are store back in data frame A.

library(geosphere)

A = data.frame(longitude = c(-2.3954998, -2.0650243, -2.0650542), latitude = c(55.32043, 55.59062, 55.60859))
B = data.frame(longitude = c(-2.4252843, -2.0650542, -2.0650243), latitude = c(55.15858, 55.60859, 55.59062))

for(i in 1:nrow(A)){
  #calucate distance against all of B
  distances<-geosphere::distGeo(A[i,], B)/1000
  #rank the calculated distances
  ranking<-rank(distances, ties.method = "first")

  #find the 3 shortest and store the indexes of B back in A
  A$shortest[i]<-which(ranking ==1) #Same as which.min()
  A$shorter[i]<-which(ranking==2)
  A$short[i]<-which(ranking ==3)

  #store the distances back in A
  A$shortestD[i]<-distances[A$shortest[i]] #Same as min()
  A$shorterD[i]<-distances[A$shorter[i]]
  A$shortD[i]<-distances[A$short[i]]
}
A

  longitude latitude shortest shorter short shortestD  shorterD   shortD
1 -2.395500 55.32043        1       3     2  18.11777 36.633310 38.28952
2 -2.065024 55.59062        3       2     1   0.00000  2.000682 53.24607
3 -2.065054 55.60859        2       3     1   0.00000  2.000682 55.05710

As M Viking pointed out, for the geosphere package the data must be arranged Lon then Lat.

Sierrasiesser

actually one last question: would you know how to display references instead of indexes in the columns shortest, shorter, short? – Apps

The index number returned is just the row number from B. You can simply take the provided indexes and reference them back to the original data frame with something like B$longitude[index] or any other column of interest. Hopefully this provided the tools you need to get to your end result with your actual data. – Sierrasiesser

Answer 2

I add below a solution using the spatialrisk package. The key functions in this package are written in C++ (Rcpp), and are therefore very fast.

The function spatialrisk::points_in_circle calculates the observations within radius from a center point. Note that distances are calculated using the Haversine formula. Since each element of the output is a data frame, purrr::map_dfr is used to row-bind them together:

purrr::map2_dfr(A$latitude, A$longitude, 
                  ~spatialrisk::points_in_circle(B, .y, .x, 
                                                 lon = longitude, 
                                                 lat = latitude, 
                                                 radius = 1e6)[1:3,], 
                .id = "id_A")

  id_A reference latitude longitude distance_m
1    1         C 55.15858 -2.425284  18115.958
2    1         E 55.59062 -2.065024  36603.447
3    1         D 55.60859 -2.065054  38260.562
4    2         E 55.59062 -2.065024      0.000
5    2         D 55.60859 -2.065054   2000.412
6    2         C 55.15858 -2.425284  53219.597
7    3         D 55.60859 -2.065054      0.000
8    3         E 55.59062 -2.065024   2000.412
9    3         C 55.15858 -2.425284  55031.092

Unscathed

much appreciate the fast code. was working on big data and the for loop one took ages. +1 – Dismantle

Answer 3

geosphere library has several functions to help you. distGeo returns meters.

Note the data must be arranged Lon then Lat.

library(geosphere)

A = data.frame(longitude = c(-2.3954998, -2.0650243, -2.0650542), latitude = c(55.32043, 55.59062, 55.60859))

B = data.frame(longitude = c(-2.4252843, -2.0650542, -2.0650243), latitude = c(55.15858, 55.60859, 55.59062))

geosphere::distGeo(A, B)

# > geosphere::distGeo(A, B)
# [1] 18117.765  2000.682  2000.682

Vector of distances in meters

Schapira

Answer 4

I know this is a long way but, in this question, there exists a formula for calculation the distance on your own. So if we convert those codes into the R we can do the same by just using base R.

Function :

rad = function(x) {
    return(x * pi / 180)

}   

getDistance = function(p1, p2) {

        R = 6378137 #  Earth’s mean radius in meter
        dLat = rad(p2[1] - p1[1])
        dLong = rad(p2[2] - p1[2])


        a = ( sin(dLat / 2) * sin(dLat / 2) +
        cos(rad(p1[1])) * cos(rad(p2[1])) *
            sin(dLong / 2) * sin(dLong / 2)  )


        c = 2 * atan2(sqrt(a),sqrt(1 - a))
        d = R * c
  return(d)  # returns the distance in meter
}

Example :

p1 <- c(55.32043 , -2.395500)
p3 <- c(55.15858 , -2.425284)

getDistance(p1,p3)
18115.96

Thus, once we can call those two functions, we can calculate any distance between two locations. So,

output <-lapply( 1:nrow(A), function(i) 
         lapply(1:nrow(B), function(j) 
             cbind(A[i,],B[j,],Distance=getDistance(as.numeric(A[i,-1]),as.numeric(B[j,-1])))

           ))

do.call(rbind,lapply(1:3,function(i) do.call(rbind,output[[i]])))

gives,

   reference latitude longitude reference latitude longitude  Distance
1          C 55.32043 -2.395500         C 55.15858 -2.425284 18115.958
2          C 55.32043 -2.395500         D 55.60859 -2.065054 38260.562
3          C 55.32043 -2.395500         E 55.59062 -2.065024 36603.447
23         D 55.59062 -2.065024         C 55.15858 -2.425284 53219.597
21         D 55.59062 -2.065024         D 55.60859 -2.065054  2000.412
22         D 55.59062 -2.065024         E 55.59062 -2.065024     0.000
33         E 55.60859 -2.065054         C 55.15858 -2.425284 55031.092
31         E 55.60859 -2.065054         D 55.60859 -2.065054     0.000
32         E 55.60859 -2.065054         E 55.59062 -2.065024  2000.412

Fai