对于表中的每个观测值，根据纬度和经度(R)计算表中x米范围内的其他观测值的数量

Question

我有一个位置及其经度和纬度的数据帧，看起来像：

set.seed(211)
latitude <-runif(5000, min=50, max=55)
longitude <- runif(5000, min=-2, max=0)
location_id <- seq(1,5000)
reprex <- data.frame(location_id, latitude, longitude)

对于每个location_id，我需要计算列表中距离该位置10英里(约16000米(以内的其他位置的数量。

我曾想过在某种for循环(或者应用函数(中使用geosphere:：distGeo((来实现这一点，但我无法解决如何对其进行编码，使其将列表中的其他项目与当前项目进行比较，并计算有多少项目在某个阈值内，记录该值，然后移动到下一行。

有人知道怎么写吗？

Answer 1

distGeo函数可以做到这一点，但您需要一个循环。请注意，坐标的第一列必须是经度。

lst <- vector(50, mode="list")
for(i in 1:50) {
    dist <- distGeo(p1=reprex[i,c(3,2)], p2=reprex[-i,c(3,2)])
    lst[[i]] <- sum(dist<16000)
}
reprex$n <- unlist(lst)
table(unlist(lst))
 0  1  2 
34 10  6

因此，50个点中有34个在10英里(约16000米(内没有其他点，10个只有1个，6个有2个。

Answer 2

fields中的rdist.earth函数似乎对此很有用，例如：

library(fields)
dist.matrix <- rdist.earth(reprex[-1])
colSums(dist.matrix<10)

在这种情况下，唯一的技巧是在逻辑矩阵上使用colSums来计数TRUE值的数量。

请注意，英里是默认值，km可以与参数miles=FALSE一起使用。

Answer 3

将循环隐藏在(仍然很慢的(apply中，并解开纬度和经度(它们通常是相反的(，您可以尝试类似的方法

set.seed(211)
latitude <-runif(5000, min=50, max=55)
longitude <- runif(5000, min=-2, max=0)
location_id <- seq(1, 5000)
reprex <- data.frame(location_id, latitude, longitude)
library(geosphere)
within10m <- function(p1, p2, dist=16000){
  sum(geosphere::distGeo(p1, p2) <= dist)
  }
matpoints <- as.matrix(reprex[, 3:2])
reprex$neighbours <- 
  apply(matpoints, 1, within10m, p2=matpoints) - 1
head(reprex) 
#   location_id latitude  longitude neighbours
# 1           1 51.17399 -1.1489713         48
# 2           2 54.52623 -1.8554624         39
# 3           3 54.84852 -0.3014742         56
# 4           4 51.72104 -1.8644226         50
# 5           5 51.32793 -0.7417923         56
# 6           6 50.07346 -0.8939857         36

Answer 4

最终我在这里使用了答案，因为它非常优雅，避免了循环：计算特定半径内的点数

我使用了代码：

library(geosphere) # for distHaversine() and distm() functions
reprex <- cbind(reprex, # appends to the dataset... 
                     count_nearby=rowSums( # ... a column which counts the rows in the dataset...
                       distm (reprex[,3:2], fun = distHaversine) # ... where the distance between current and other rows...
                       <= 16000)-1 # ... is less than 16000 metres. Take one away because it counts itself!
                ) #close the cbind brackets!