lat/lon 点（坐标）的 2 个列表之间的地理/地理空间距离
Geographic / geospatial distance between 2 lists of lat/lon points (coordinates)
我有 2 个列表（list1、list2），其中包含各个位置的纬度/经度。一个列表 (list2) 包含 list1 没有的地方名称。
我也想要 list1 中每个点的大概位置。所以我想在 list1 中取一个点，尝试在 list2 中寻找最近的点并取那个地方。我重复 list1 中的每个点。它还需要距离（以米为单位）和点的索引（以 list1 为单位），因此我可以围绕它构建一些业务规则 - 本质上，这些是应该添加到 list1 的 2 个新列（ 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"))
要计算具有 latitude/longitude 坐标的两点之间的地理距离，您可以使用多个公式。 geosphere 包中有 distCosine、distHaversine、distVincentySphere 和 distVincentyEllipsoid 用于计算距离。其中，distVincentyEllipsoid 被认为是最准确的，但计算量比其他的要大。
使用这些函数之一，您可以制作距离矩阵。基于该矩阵，您可以根据 which.min 的最短距离和 min 的相应距离分配 locality 名称（参见答案的最后一部分），如下所示：
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
如您所见，这在大多数情况下（十分之七）导致另一个分配 locality。
您可以添加距离：
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
感谢 Martin Haringa 提供的此解决方案，当您需要通过遍历 Mark Needham's blog
上的数据框来执行此功能时，该解决方案使这种方式变得更容易
library(dplyr)
library(geosphere)
df %>%
rowwise() %>%
mutate(newcolumn_distance = distHaversine(c(df$long1, df$lat1),
c(df$long2, df$lat2)))
我在真实世界数据集的大样本上分别使用 distm 和 distHaversine 这两个函数进行了测试，distHaversine 似乎比 distm 函数快得多。我很惊讶，因为我认为这两者只是两种格式的相同功能。
我在下面添加了一个使用 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
我有 2 个列表（list1、list2），其中包含各个位置的纬度/经度。一个列表 (list2) 包含 list1 没有的地方名称。
我也想要 list1 中每个点的大概位置。所以我想在 list1 中取一个点，尝试在 list2 中寻找最近的点并取那个地方。我重复 list1 中的每个点。它还需要距离（以米为单位）和点的索引（以 list1 为单位），因此我可以围绕它构建一些业务规则 - 本质上，这些是应该添加到 list1 的 2 个新列（ 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"))
要计算具有 latitude/longitude 坐标的两点之间的地理距离，您可以使用多个公式。 geosphere 包中有 distCosine、distHaversine、distVincentySphere 和 distVincentyEllipsoid 用于计算距离。其中，distVincentyEllipsoid 被认为是最准确的，但计算量比其他的要大。
使用这些函数之一，您可以制作距离矩阵。基于该矩阵，您可以根据 which.min 的最短距离和 min 的相应距离分配 locality 名称（参见答案的最后一部分），如下所示：
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
如您所见，这在大多数情况下（十分之七）导致另一个分配 locality。
您可以添加距离：
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
感谢 Martin Haringa 提供的此解决方案，当您需要通过遍历 Mark Needham's blog
上的数据框来执行此功能时，该解决方案使这种方式变得更容易library(dplyr)
library(geosphere)
df %>%
rowwise() %>%
mutate(newcolumn_distance = distHaversine(c(df$long1, df$lat1),
c(df$long2, df$lat2)))
我在真实世界数据集的大样本上分别使用 distm 和 distHaversine 这两个函数进行了测试，distHaversine 似乎比 distm 函数快得多。我很惊讶，因为我认为这两者只是两种格式的相同功能。
我在下面添加了一个使用 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