两个经纬度点列表之间的地理/地理空间距离

Question

我有两个列表（list1，list2），分别包含不同位置的纬度/经度坐标。其中一个列表（list2）包含了list1没有的地名。

我想为list1中的每个点找到大致的所在地。因此，我想要取list1中的一个点，尝试查找list2中最近的点，并获取该地名。我对list1中的每个点都要重复此操作。我还需要距离（以米为单位）和该点在list1中的索引，以便我可以建立一些业务规则 - 本质上这是应添加到list1中的两个新列 (near_dist, indx)。

我正在使用gdist函数，但我无法将其与数据框输入一起使用。

示例输入列表：

list1 <- data.frame(longitude = c(80.15998, 72.89125, 77.65032, 77.60599, 
                                  72.88120, 76.65460, 72.88232, 77.49186, 
                                  72.82228, 72.88871), 
                    latitude = c(12.90524, 19.08120, 12.97238, 12.90927, 
                                 19.08225, 12.81447, 19.08241, 13.00984,
                                 18.99347, 19.07990))
list2 <- data.frame(longitude = c(72.89537, 77.65094, 73.95325, 72.96746, 
                                  77.65058, 77.66715, 77.64214, 77.58415,
                                  77.76180, 76.65460), 
                    latitude = c(19.07726, 13.03902, 18.50330, 19.16764, 
                                 12.90871, 13.01693, 13.00954, 12.92079,
                                 13.02212, 12.81447), 
                    locality = c("A", "A", "B", "B", "C", "C", "C", "D", "D", "E"))

- grammar

Answer 1

要计算具有纬度/经度坐标的两点之间的地理距离，可以使用几个公式。包geosphere具有distCosine、distHaversine、distVincentySphere和distVincentyEllipsoid等用于计算距离的函数。其中，distVincentyEllipsoid被认为是最准确的，但比其他函数更需要计算。

使用其中一个函数，您可以创建一个距离矩阵。然后可以基于该矩阵使用which.min指定locality名称的最短距离，并使用min指定相应的距离（请参见答案的最后一部分），如下所示：

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)]

这是什么意思：

> 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

---

另一种可能的方法是根据list2中locality的平均经纬度值来分配locality:

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)])

或者使用 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)] ]

这将会得到：

> 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

如您所见，这在大多数情况下（10次中有7次）会导致另一个指定的地区。

---

您可以通过以下方式添加距离：

list1$near_dist <- apply(mat2, 1, min)

或者使用 max.col 的另一种方法（这种方法很可能更快）：

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)] ]

结果：

> 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

- Jaap

Answer 2

感谢Martin Haringa提供的解决方案，在Mark Needham的博客上学习如何更轻松地通过遍历数据框执行此函数。

library(dplyr)
library(geosphere)

df %>%
  rowwise() %>%
  mutate(newcolumn_distance = distHaversine(c(df$long1, df$lat1), 
                                            c(df$long2, df$lat2)))

我在真实世界数据集的大样本上分别使用了distm和distHaversine两个函数进行测试，结果distHaversine似乎比distm函数快得多。我感到惊讶，因为我认为这两个函数只是以两种不同的格式表示相同的函数。

- Alexander Kielland

---

博客链接已损坏。 - Jeff Parker

抱歉 @JeffParker，现在应该已经修复了！ - Alexander Kielland

@AlexanderKielland 谢谢。更简单的方法是使用矢量化版本的haversine函数，例如：df％>% mutate（newcolumn_distance = spatialrisk :: haversine（lat1，long1，lat2，long2）） - mharinga

不错的@mharinga！我有机会时一定要测试一下！我有一个巨大的数据集，肯定需要矢量化选项来提高吞吐量！ - Alexander Kielland

Answer 3

我使用spatialrisk包提供以下解决方案。该包中的关键函数是用C++（Rcpp）编写的，因此速度非常快。

函数spatialrisk::points_in_circle()计算距离中心点半径范围内的观测值。请注意，距离是使用Haversine公式计算的。由于输出的每个元素都是数据框，因此使用purrr::map_dfr将它们行绑在一起：

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

- mharinga