Для каждого наблюдения в таблице подсчитайте количество других наблюдений в таблице в пределах x метров на основе широты и долготы (R)

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 метров) этого места.

Я думал об использовании geosphere :: distGeo () для этого в некотором роде для l oop (или, возможно, функции применения), но я просто не могу работать как его кодировать, чтобы он сравнивал каждый второй элемент в списке с текущим элементом и подсчитывал, сколько их находится в пределах определенного порога, записывает это значение и затем переходит к следующей строке.

Кто-нибудь знает как это написать?

Answer 1

В конечном итоге я использовал ответ здесь, так как он был довольно элегантным и избегал циклов: Вычисление количества точек в пределах определенного радиуса

Я использовал код:

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!

Answer 2

Функция rdist.earth в fields кажется полезной для этого, например:

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

Единственный прием в этом случае - использование colSums в логической матрице для подсчета числа TRUE значения.

Обратите внимание, что по умолчанию используются мили, км можно использовать с аргументом miles=FALSE.

Answer 3

Скрывая l oop в (все еще медленном) 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

Функция distGeo может сделать это, но вам нужно al oop. Обратите внимание, что первый столбец координат должен быть долготой.

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

Таким образом, 34 из 50 точек не имеют никакой другой точки в пределах 10 миль (~ 16 000 метров), 10 имеют только 1, а 6 имеют 2.