5

我有两个数据集 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 中这些点的经纬度的参考和坐标。

我尝试使用这篇文章

R - 查找给定半径内最近的相邻点和相邻点的数量，坐标经纬度

但是，即使我按照所有说明进行操作，主要是使用distmpackage中的命令geosphere，距离仍以不可能为公里的单位出现。我看不出代码中有什么要改变的，特别是因为我对这些geo包一点也不熟悉。

4

4 回答 4

2

我在下面添加了一个使用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
于 2019-10-25T14:19:53.760 回答
1

geosphere图书馆有几个功能可以帮助你。distGeo返回米。

注意数据必须排列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

以米为单位的距离向量

于 2019-08-16T14:10:14.467 回答
1

这是使用单个循环和矢量化距离计算（转换为公里）的解决方案。
该代码使用base R的rank函数对计算距离列表进行排序/排序。
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 包，数据必须按 Lon 然后 Lat 排列。

于 2019-08-16T14:45:45.430 回答
0

我知道这是一条很长的路，但是在这个问题中，有一个公式可以自己计算距离。因此，如果我们将这些代码转换为 ，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
于 2019-08-16T14:33:19.460 回答