r - 识别R中指定距离内的点

Question

从包含十进制坐标的 10,000 个站点的列表中，我尝试根据这些站点之间计算的距离来识别彼此相距在 100 英尺以内的站点，并创建这些站点的子集。在最终列表中，我希望获得彼此相距 100 英尺以内的车站名称、它们的纬度和经度以及它们之间的距离。

我在 mathworks(使用 rangesearch)或 SQL 或 JAVA 等其他平台上发现了类似的问题，但在 R 中没有发现类似的问题。

R 有办法做到这一点吗？我找到的最接近的答案是 Listing number of obervations by location其中列出了一定距离内的观测数量，但似乎答案不完整，无法确定彼此之间特定距离内的站点。

基本上我想找出哪些站位于同一地点。

我非常感谢任何对此的帮助。

Answer 1

两种方法。

第一个使用fossil包中的earth.dist(...)创建一个距离矩阵，然后利用data.tables来组装表格结果。

第二个方法使用 geosphere 包中的 distHaversine(...) 来计算距离并一步组装最终的共置表。后一种方法可能会更快，也可能不会更快，但肯定会更有效地存储内存，因为它从不存储完整的距离矩阵。此外，这种方法还可以在地理圈中使用其他距离测量，例如，distVincentySphere(...)、distVincentyEllipsoid(...) ，或distMeeus(...)。

请注意，实际距离略有不同，可能是因为 earth.dist(...) 和 distHaversine(...) 使用的半径估计值略有不同地球的。另请注意，这里的两种方法都依赖于 ID 的站号。如果电台有名称，则代码需要稍作修改。

第一种方法:使用earth.dist(...)

df = read.table(header=T,text="long lat
                1 -74.20139 39.82806
                2 -74.20194 39.82806 
                3 -74.20167 39.82806 
                4 -74.20197 39.82824 
                5 -74.20150 39.82814 
                6 -74.26472 39.66639 
                7 -74.17389 39.87111 
                8 -74.07224 39.97353 
                9 -74.07978 39.94554")              # your sample data
library(fossil)                                     # for earth.dist(...)
library(data.table)
sep.ft   <- 200                                     # critical separation (feet)
sep.km   <- sep.ft*0.0003048                        # critical separation (km)
m        <- as.matrix(earth.dist(df))               # distance matrix in km
coloc    <- data.table(which(m<sep.km, arr.ind=T))  # pairs of stations with dist<200 ft
setnames(coloc,c("row","col"),c("ST.1","ST.2"))     # rename columns to reflect station IDs
coloc    <- coloc[ST.1<ST.2,]                       # want only lower triagular part
coloc[,dist:=m[ST.1,ST.2]/0.0003048,by="ST.1,ST.2"] # append distances in feet
remove(m)                                           # don't need distance matrix anymore...
stations <- data.table(id=as.integer(rownames(df)),df)
setkey(stations,id)
setkey(coloc,ST.1)
coloc[stations,c("long.1","lat.1"):=list(long,lat),nomatch=0]
setkey(coloc,ST.2)
coloc[stations,c("long.2","lat.2"):=list(long,lat),nomatch=0]

产生这个:

coloc
#     ST.1 ST.2      dist    long.1    lat.1    long.2    lat.2
#  1:    1    2 154.13436 -74.20139 39.82806 -74.20194 39.82806
#  2:    1    3  78.46840 -74.20139 39.82806 -74.20167 39.82806
#  3:    2    3  75.66596 -74.20194 39.82806 -74.20167 39.82806
#  4:    1    4 175.31180 -74.20139 39.82806 -74.20197 39.82824
#  5:    2    4  66.22069 -74.20194 39.82806 -74.20197 39.82824
#  6:    3    4 106.69018 -74.20167 39.82806 -74.20197 39.82824
#  7:    1    5  42.45634 -74.20139 39.82806 -74.20150 39.82814
#  8:    2    5 126.71608 -74.20194 39.82806 -74.20150 39.82814
#  9:    3    5  55.87449 -74.20167 39.82806 -74.20150 39.82814
# 10:    4    5 136.67612 -74.20197 39.82824 -74.20150 39.82814

第二种方法:使用distHaversine(...)

library(data.table)
library(geosphere)
sep.ft   <- 200                       # critical separation (feet)
stations <- data.table(id=as.integer(rownames(df)),df)

d <- function(x){                     # distance between station[i] and all subsequent stations
  r.ft <- 6378137*3.28084             # radius of the earth, in feet
  if (x[1]==nrow(stations)) return()  # don't process last row
  ref <- stations[(x[1]+1):nrow(stations),]
  z <- distHaversine(ref[,2:3,with=F],x[2:3], r=r.ft)
  z <- data.table(ST.1=x[1], ST.2=ref$id, dist=z, long.1=x[2], lat.1=x[3], long.2=ref$long, lat.2=ref$lat)
  return(z[z$dist<sep.ft,])
}
coloc.2 = do.call(rbind,apply(stations,1,d))

产生这个:

coloc.2
#     ST.1 ST.2      dist    long.1    lat.1    long.2    lat.2
#  1:    1    2 154.26350 -74.20139 39.82806 -74.20194 39.82806
#  2:    1    3  78.53414 -74.20139 39.82806 -74.20167 39.82806
#  3:    1    4 175.45868 -74.20139 39.82806 -74.20197 39.82824
#  4:    1    5  42.49191 -74.20139 39.82806 -74.20150 39.82814
#  5:    2    3  75.72935 -74.20194 39.82806 -74.20167 39.82806
#  6:    2    4  66.27617 -74.20194 39.82806 -74.20197 39.82824
#  7:    2    5 126.82225 -74.20194 39.82806 -74.20150 39.82814
#  8:    3    4 106.77957 -74.20167 39.82806 -74.20197 39.82824
#  9:    3    5  55.92131 -74.20167 39.82806 -74.20150 39.82814
# 10:    4    5 136.79063 -74.20197 39.82824 -74.20150 39.82814