Работа с координатами и огромными наборами данных в R

Question

Я немного борюсь с двумя наборами данных, содержащими координаты людей и вышек сотовой связи:

• Первый набор данных по 9 459 лицам с 1214 переменными, включая их широту и долготу в градусах.
• второй набор данных о 31 176 сотовых вышках с 4 переменными, включая их широту и долготу в градусах и диапазон в метрах.

Я хотел бы определить, находится ли человек в диапазоне хотя бы одной из вышек сотовой связи, и создать пустышку, равную 1, если это так.

Однако из-за размера наборов данных я не могу объединить их с помощью команды перекрестного соединения. Я пытался использовать geosphere пакет с помощью следующей команды:

distm(c(df1$longitude, df2$latitude), c(df2$longitude, df2$latitude), fun= distHaversine)

К сожалению, это не работает, так как два набора данных не имеют одинакового размера. Есть идеи, как решить эту проблему?

Answer 1

Ниже я добавляю решение с использованием пакета пространственного риска. Ключевые функции в этом пакете написаны на C++ (Rcpp) и поэтому работают очень быстро.

Функция Spacerisk::points_in_circle() вычисляет наблюдения в радиусе от центральной точки. Обратите внимание, что расстояния рассчитываются по формуле Хаверсина. Поскольку каждый элемент вывода - это фрейм данных, purrr::map_dfr используется для их строкового связывания:

library(tibble)
library(spatialrisk)
library(dplyr)

set.seed(1702)
users <- tibble(userid = as.character(1:10000),
                lon = rnorm(10000, 16.3738, 1),
                lat = rnorm(10000, 48.2082, 1))

towers <- tibble(lon = rnorm(35000, 16.3738, 1),
                 lat = rnorm(35000, 48.2082, 1))

# Users with tower within 200 meters
purrr::map2_dfr(users$lon, users$lat, 
                   ~points_in_circle(towers, .x, .y, radius = 200)[1,], 
                   .id = "userid") %>%
     mutate(inrange = ifelse(is.na(distance_m), FALSE, TRUE))

Answer 2

Как правило, это можно сделать гораздо эффективнее, чтобы максимально использовать оперативную память и процессор и снизить накладные расходы. Однако, если то, что вы пытаетесь сделать, является одноразовой операцией, подходящего ниже подхода должно быть достаточно (на текущем ноутбуке требуется около 5 минут).

Вспомогательная функция

# More info: https://github.com/RomanAbashin/distGeo_v
distGeo_v <- function(x, y, xx, yy) { 
    if(!"geosphere" %in% installed.packages())  {
        stop("The 'geosphere' package needs to be installed for this function to work.")
    }
    matrix(.Call("_inversegeodesic", 
                 as.double(x), as.double(y), as.double(xx), as.double(yy), 
                 as.double(6378137), 1/298.257223563, PACKAGE='geosphere'), 
           ncol = 3, byrow = TRUE)[,1]
}

Данные

library(geosphere)
library(tidyverse)
set.seed(1702)

users <- tibble(userid = 1:10000,
                x = rnorm(10000, 16.3738, 5),
                y = rnorm(10000, 48.2082, 5))

towers <- tibble(lon = rnorm(35000, 16.3738, 10),
                 lat = rnorm(35000, 48.2082, 10),
                 range = runif(35000, 50, 10000))

Код

result <- NULL
for(i in 1:nrow(users)) {

    is_match <- users[i, 1:3] %>%
        tidyr::crossing(towers[, 1:3]) %>%
        filter(distGeo_v(x, y, lon, lat) <= range) %>%
        nrow() > 0

    result <- bind_rows(result, tibble(userid = users$userid[i],
                                       match = is_match))

}

Результат

> head(result)
# A tibble: 6 x 2
  userid match
   <int> <lgl>
1      1 TRUE 
2      2 FALSE
3      3 FALSE
4      4 TRUE 
5      5 FALSE
6      6 FALSE

Теперь вы можете left_join результат к вашим исходным данным.