检查坐标是否在给定半径内

Question

我有一个这种格式的某些巴士站的坐标列表

 Bus_Stop_ID     lat          long
 A               -34.04199    18.61747
 B               -33.92312    18.44649

然后我有一些商店的列表

 Shop_ID     lat          long
 1            -34.039350  18.617964  
 2            -33.927820  18.410520 

我想检查商店是否在巴士站方圆 500 米范围内。最终，最终数据集将如下所示，其中 Bus_Stop 列指示 T/F，如果 Bus_Stop == T -

，则 Bus_Stop_ID 显示该商店的相关 BUS ID

 Shop_ID     lat          long       Bus_Stop Bus_ID  
 1            -34.039350  18.617964  TRUE     A
 2            -33.927820  18.410520  FALSE    #NA

有人知道我如何使用 R 来解决这个问题吗？我已经看过 geosphere 包，但由于我在空间领域的相对缺乏经验，我很难理解它。您可以推荐任何想法或软件包吗？谢谢

Answer 1

更新为更具可扩展性的解决方案：

前面的答案（仍然包含在下面）不适合大型数据集。原因是我们需要计算每对shops 和bus 的距离。因此，N 商店和M 总线的内存和计算规模都为O(N*M)。更具可扩展性的解决方案使用 KD-Tree 等数据结构对每个商店执行最近邻搜索。这里的优点是计算复杂度变为O(M*logM) 用于构建公交车站的KD-Tree 和O(N*logM) 用于搜索每个商店的最近邻居。

为此，我们可以使用 RANN 包中的 nn2。这里的复杂之处在于nn2 仅处理欧几里得距离，并且不知道任何关于纬度/经度的信息。因此，我们需要将纬度/经度坐标转换为一些地图投影（即 UTM），以便正确使用它（即，为了正确计算商店和公交车站之间的欧几里得距离）。

注意：以下内容大量借鉴了 Josh O'Brien 针对 determining the UTM zone from a longitude 和 converting lat/long to UTM 的解决方案，因此他应该鞠躬。

## First define a function from Josh OBrien's answer to convert
## a longitude to its UTM zone
long2UTM <- function(long) {
  (floor((long + 180)/6) %% 60) + 1
}

## Assuming that all points are within a zone (within 6 degrees in longitude),
## we use the first shop's longitude to get the zone.
z <- long2UTM(shops[1,"long"])

library(sp)
library(rgdal)

## convert the bus lat/long coordinates to UTM for the computed zone
## using the other Josh O'Brien linked answer
bus2 <- bus
coordinates(bus2) <- c("long", "lat")
proj4string(bus2) <- CRS("+proj=longlat +datum=WGS84")

bus.xy <- spTransform(bus2, CRS(paste0("+proj=utm +zone=",z," ellps=WGS84")))

## convert the shops lat/long coordinates to UTM for the computed zone
shops2 <- shops
coordinates(shops2) <- c("long", "lat")
proj4string(shops2) <- CRS("+proj=longlat +datum=WGS84")

shops.xy <- spTransform(shops2, CRS(paste0("+proj=utm +zone=",z," ellps=WGS84")))

library(RANN)

## find the nearest neighbor in bus.xy@coords for each shops.xy@coords
res <- nn2(bus.xy@coords, shops.xy@coords, 1)
## res$nn.dist is a vector of the distance to the nearest bus.xy@coords for each shops.xy@coords
## res$nn.idx is a vector of indices to bus.xy of the nearest bus.xy@coords for each shops.xy@coords
shops$Bus_Stop <- res$nn.dists <= 500
shops$Bus_ID <- ifelse(res$nn.dists <= 500, bus[res$nn.idx,"Bus_Stop_ID"], NA)

虽然更复杂，但这种方法更适合处理可能有大量商店和公交车站的实际问题。使用相同的提供数据：

print(shops)
##  Shop_ID       lat     long Bus_Stop Bus_ID
##1       1 -34.03935 18.61796     TRUE      A
##2       2 -33.92782 18.41052    FALSE   <NA>

---

您可以使用包geosphere 执行此操作。在这里，我假设您的第一个数据框命名为bus，而您的第二个数据框命名为shops：

library(geosphere)
g <- expand.grid(1:nrow(shops), 1:nrow(bus))
d <- matrix(distGeo(shops[g[,1],c("long","lat")], bus[g[,2],c("long","lat")]),
            nrow=nrow(shops))
shops$Bus_Stop <- apply(d, 1, function(x) any(x <= 500))
shops$Bus_ID <- bus[apply(d, 1, function(x) {
                                  c <-which(x <= 500)
                                  if(length(c)==0) NA else c[1]
                                }), "Bus_Stop_ID"]
print(shops)
##  Shop_ID       lat     long Bus_Stop Bus_ID
##1       1 -34.03935 18.61796     TRUE      A
##2       2 -33.92782 18.41052    FALSE   <NA>

注意事项：

1. 我们首先使用expand.grid 来枚举shops 和bus 停靠点的所有对组合。这些是由shops 首先订购的。
2. 然后我们使用geosphere::distGeo 计算距离矩阵d。请注意，输入需要 (lon, lat) 坐标。 distGeo 返回以米为单位的距离。生成的 d 矩阵是 now(shops) 乘以 now(bus)，因此每一行给出了从商店到每个公共汽车站的距离。
3. 然后我们通过在d 中使用apply 和MARGIN=1 对每一行x 应用函数any(x <= 500) 来查看每个商店500 米内是否有公交车站。
4. 同样，我们可以在我们应用的函数中使用which而不是any提取500米内的第一家店铺d的列（对应bus中的行）。然后使用此结果从bus 中选择Bus_Stop_ID。

顺便说一句，我们不必apply 条件x <= 500 两次。以下也将起作用：

shops$Bus_ID <- bus[apply(d, 1, function(x) {
                                  c <-which(x <= 500)
                                  if(length(c)==0) NA else c[1]
                                }), "Bus_Stop_ID"]
shops$Bus_Stop <- !is.na(shops$Bus_ID)

而且效率更高。

数据：

bus <- structure(list(Bus_Stop_ID = structure(1:2, .Label = c("A", "B"
), class = "factor"), lat = c(-34.04199, -33.92312), long = c(18.61747, 
18.44649)), .Names = c("Bus_Stop_ID", "lat", "long"), class = "data.frame",  row.names = c(NA, 
-2L))

shops <- structure(list(Shop_ID = 1:2, lat = c(-34.03935, -33.92782), 
long = c(18.617964, 18.41052), Bus_ID = structure(c(1L, NA
), .Label = c("A", "B"), class = "factor"), Bus_Stop = c(TRUE, 
FALSE)), .Names = c("Shop_ID", "lat", "long", "Bus_ID", "Bus_Stop"
), row.names = c(NA, -2L), class = "data.frame")

Answer 2

我的第一种方法是只使用Euclidean distance 并检查结果值是否大于或等于 0。

然后您可以使用 IF 子句并检查 T/F 条件。

我希望这会有所帮助。

PS：在我的想象中，500m 的距离将是地球表面相当平坦的表示，所以我认为不需要使用一些大地水准面包。