查找从一个数据集到第二个数据集的最近点(纬度/经度)
Find closest points (lat / lon) from one data set to a second data set
我有两个数据集,A 和 B,它们给出了英国不同点的位置:
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 有 400 行,B 有 1800 行。
对于 A 中的所有行,我想找到 A 中的一个点与 B 中三个最近点中的每一个之间的最短距离(以千米为单位),以及 B 中这些点的经纬度参考和坐标。
我试过用这个post
R - Finding closest neighboring point and number of neighbors within a given radius, coordinates lat-long
然而,即使我按照所有说明进行操作,主要使用包 geosphere
中的命令 distm
,距离的单位也不可能是公里。我看不到要更改代码的内容,尤其是因为我对 geo
包一点都不熟悉。
geosphere
图书馆有几个功能可以帮助您。 distGeo
returns米。
注意数据必须排列Lon
然后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
以米为单位的距离向量
我知道这条路很远,但是在this问题中,有一个公式可以自己计算距离。因此,如果我们将这些代码转换为 R
,我们只需使用 base R
即可完成相同的操作。
函数:
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
}
示例:
p1 <- c(55.32043 , -2.395500)
p3 <- c(55.15858 , -2.425284)
getDistance(p1,p3)
18115.96
因此,一旦我们可以调用这两个函数,我们就可以计算两个位置之间的任何距离。所以,
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]])))
给予,
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
这是使用单个循环并矢量化距离计算(转换为公里)的解决方案。
该代码使用基础 R 的 rank
函数来 order/sort 计算距离列表。
索引和计算出的 3 个最短值的距离存储回数据框 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
正如 M Viking 指出的那样,对于 geosphere 包,数据必须先经纬度排列。
我在下面添加了一个使用 spatialrisk
包的解决方案。这个包中的关键函数是用 C++ (Rcpp) 编写的,因此速度非常快。
函数spatialrisk::points_in_circle
计算中心点半径范围内的观测值。请注意,距离是使用 Haversine 公式计算的。由于输出的每个元素都是一个数据框,因此 purrr::map_dfr
用于将它们行绑定在一起:
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
我有两个数据集,A 和 B,它们给出了英国不同点的位置:
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 有 400 行,B 有 1800 行。
对于 A 中的所有行,我想找到 A 中的一个点与 B 中三个最近点中的每一个之间的最短距离(以千米为单位),以及 B 中这些点的经纬度参考和坐标。
我试过用这个post
R - Finding closest neighboring point and number of neighbors within a given radius, coordinates lat-long
然而,即使我按照所有说明进行操作,主要使用包 geosphere
中的命令 distm
,距离的单位也不可能是公里。我看不到要更改代码的内容,尤其是因为我对 geo
包一点都不熟悉。
geosphere
图书馆有几个功能可以帮助您。 distGeo
returns米。
注意数据必须排列Lon
然后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
以米为单位的距离向量
我知道这条路很远,但是在this问题中,有一个公式可以自己计算距离。因此,如果我们将这些代码转换为 R
,我们只需使用 base R
即可完成相同的操作。
函数:
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
}
示例:
p1 <- c(55.32043 , -2.395500)
p3 <- c(55.15858 , -2.425284)
getDistance(p1,p3)
18115.96
因此,一旦我们可以调用这两个函数,我们就可以计算两个位置之间的任何距离。所以,
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]])))
给予,
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
这是使用单个循环并矢量化距离计算(转换为公里)的解决方案。
该代码使用基础 R 的 rank
函数来 order/sort 计算距离列表。
索引和计算出的 3 个最短值的距离存储回数据框 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
正如 M Viking 指出的那样,对于 geosphere 包,数据必须先经纬度排列。
我在下面添加了一个使用 spatialrisk
包的解决方案。这个包中的关键函数是用 C++ (Rcpp) 编写的,因此速度非常快。
函数spatialrisk::points_in_circle
计算中心点半径范围内的观测值。请注意,距离是使用 Haversine 公式计算的。由于输出的每个元素都是一个数据框,因此 purrr::map_dfr
用于将它们行绑定在一起:
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