Geographic / geospatial distance between 2 lists of lat/lon points (coordinates)
Credits to Martin Haringa for this solution on making this way easier when you need this function performed by traversing down a data frame on Mark Needham's blog
library(dplyr)
library(geosphere)
df %>%
rowwise() %>%
mutate(newcolumn_distance = distHaversine(c(df$long1, df$lat1),
c(df$long2, df$lat2)))
I tested using the two functions distm and distHaversine separately on large samples from real world datasets, and distHaversine seems to come out far faster than the distm function. I'm surprised as I thought the two were merely the same function in two formats.
To calculate the geographic distance between two points with latitude/longitude coordinates, you can use several formula's. The package geosphere has the distCosine, distHaversine, distVincentySphere and distVincentyEllipsoid for calculating the distance. Of these, the distVincentyEllipsoid is considered the most accurate one, but is computationally more intensive than the other ones.
With one of these functions, you can make a distance matrix. Based on that matrix you can then assign locality names based on shortest distance with which.min and the corresponding distance with min (see for this the last part of the answer) like this:
library(geosphere)
# create distance matrix
mat <- distm(list1[,c('longitude','latitude')], list2[,c('longitude','latitude')], fun=distVincentyEllipsoid)
# assign the name to the point in list1 based on shortest distance in the matrix
list1$locality <- list2$locality[max.col(-mat)]
this gives:
> list1 longitude latitude locality 1 80.15998 12.90524 D 2 72.89125 19.08120 A 3 77.65032 12.97238 C 4 77.60599 12.90927 D 5 72.88120 19.08225 A 6 76.65460 12.81447 E 7 72.88232 19.08241 A 8 77.49186 13.00984 D 9 72.82228 18.99347 A 10 72.88871 19.07990 A
Another possibility is to assign the locality based on the average longitude and latitude values of the localitys in list2:
library(dplyr)
list2a <- list2 %>% group_by(locality) %>% summarise_each(funs(mean)) %>% ungroup()
mat2 <- distm(list1[,c('longitude','latitude')], list2a[,c('longitude','latitude')], fun=distVincentyEllipsoid)
list1 <- list1 %>% mutate(locality2 = list2a$locality[max.col(-mat2)])
or with data.table:
library(data.table)
list2a <- setDT(list2)[,lapply(.SD, mean), by=locality]
mat2 <- distm(setDT(list1)[,.(longitude,latitude)], list2a[,.(longitude,latitude)], fun=distVincentyEllipsoid)
list1[, locality2 := list2a$locality[max.col(-mat2)] ]
this gives:
> list1 longitude latitude locality locality2 1 80.15998 12.90524 D D 2 72.89125 19.08120 A B 3 77.65032 12.97238 C C 4 77.60599 12.90927 D C 5 72.88120 19.08225 A B 6 76.65460 12.81447 E E 7 72.88232 19.08241 A B 8 77.49186 13.00984 D C 9 72.82228 18.99347 A B 10 72.88871 19.07990 A B
As you can see, this leads in most (7 out of 10) occasions to another assigned locality.
You can add the distance with:
list1$near_dist <- apply(mat2, 1, min)
or another approach with max.col (which is highly probable faster):
list1$near_dist <- mat2[matrix(c(1:10, max.col(-mat2)), ncol = 2)]
# or using dplyr
list1 <- list1 %>% mutate(near_dist = mat2[matrix(c(1:10, max.col(-mat2)), ncol = 2)])
# or using data.table (if not already a data.table, convert it with 'setDT(list1)' )
list1[, near_dist := mat2[matrix(c(1:10, max.col(-mat2)), ncol = 2)] ]
the result:
> list1 longitude latitude locality locality2 near_dist 1: 80.15998 12.90524 D D 269966.8970 2: 72.89125 19.08120 A B 65820.2047 3: 77.65032 12.97238 C C 739.1885 4: 77.60599 12.90927 D C 9209.8165 5: 72.88120 19.08225 A B 66832.7223 6: 76.65460 12.81447 E E 0.0000 7: 72.88232 19.08241 A B 66732.3127 8: 77.49186 13.00984 D C 17855.3083 9: 72.82228 18.99347 A B 69456.3382 10: 72.88871 19.07990 A B 66004.9900
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:
ans <- purrr::map2_dfr(list1$longitude,
list1$latitude,
~spatialrisk::points_in_circle(list2, .x, .y,
lon = longitude,
lat = latitude,
radius = 2000000)[1,])
cbind(list1, ans)
longitude latitude longitude latitude locality distance_m
1 80.15998 12.90524 77.76180 13.02212 D 260484.0591
2 72.89125 19.08120 72.89537 19.07726 A 616.6369
3 77.65032 12.97238 77.64214 13.00954 C 4230.7216
4 77.60599 12.90927 77.58415 12.92079 D 2694.4566
5 72.88120 19.08225 72.89537 19.07726 A 1590.8723
6 76.65460 12.81447 76.65460 12.81447 E 0.0000
7 72.88232 19.08241 72.89537 19.07726 A 1487.8028
8 77.49186 13.00984 77.58415 12.92079 D 14089.1051
9 72.82228 18.99347 72.89537 19.07726 A 12089.6454
10 72.88871 19.07990 72.89537 19.07726 A 759.8012