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

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_ID 显示该商店的相关公交车 ID(如果 Bus_Stop == T -

)

 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 。因此内存和计算规模均为 O(N*M)对于 N商店和M公共(public)汽车。更具可扩展性的解决方案使用 KD 树等数据结构来对每个商店执行最近邻居搜索。这里的优点是计算复杂度变成O(M*logM)用于为公交车站构建 KD 树和 O(N*logM)用于搜索每个商店的最近邻居。

为此，我们可以使用 nn2来自RANN包裹。这里的复杂性是nn2仅处理欧几里德距离，不知道有关纬度/经度的任何信息。因此，我们需要将纬度/经度坐标转换为某种 map 投影(即 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. 然后我们计算距离矩阵 d使用geosphere::distGeo 。请注意，此处输入需要 (lon, lat) 坐标。 distGeo返回以米为单位的距离。由此产生d矩阵是now(shops)通过now(bus)这样每一行都会给出从商店到每个公交车站的距离。
3. 然后，我们通过应用函数 any(x <= 500) 来查看每个商店 500 米范围内是否有公交车站每行 x在d使用apply与 MARGIN=1 .
4. 同样，我们可以提取d这一列(对应于 bus 中的行)使用 which 500 米内的第一家商店而不是any在我们的应用函数中。然后使用此结果选择 Bus_Stop_ID来自bus .

顺便说一句，我们不必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")