Географическое расстояние по группе — применение функции к каждой паре строк

Мои 10 центов. Ты сможешь:

# subset the province
df1 <- df1[which(df1$province==1),]

# get all combinations
all <- combn(df1$house, 2, FUN = NULL, simplify = TRUE)

# run your function and get distances for all combinations
distances <- c()
for(col in 1:ncol(all)) {
  a <- all[1, col]
  b <- all[2, col]
  dist <- distm(c(df1$lon[a], df1$lat[a]), c(df1$lon[b], df1$lat[b]), fun = distHaversine)
  distances <- c(distances, dist)
  }

# calculate mean:
mean(distances)
# [1] 15379.21

Это дает вам среднее значение для провинции, которое вы можете сравнить с результатами других методов. Например sapply который упоминался в комментариях:

df1 <- df1[which(df1$province==1),]
mean(sapply(split(df1, df1$province), dist))
# [1] 1.349036

Как видите, это дает разные результаты, потому что функция dist может вычислять расстояния разного типа (например, евклидовы) и не может вычислять гаверсинус или другие «геодезические» расстояния. В пакете geodist есть варианты, которые могут приблизить вас к sapply:

library(geodist)
library(magrittr)

# defining the data
df1 <- data.frame(province = c(1, 1, 1, 2, 2, 2),
                  house = c(1, 2, 3, 4, 5, 6),
                  lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7), 
                  lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2))

# defining the function 
give_distance <- function(resultofsplit){
  distances <- c()
  for (i in 1:length(resultofsplit)){
    sdf <- resultofsplit
    sdf <- sdf[[i]]
    sdf <- sdf[c("lon", "lat", "province", "house")]

    sdf2 <- as.matrix(sdf)
    sdf3 <- geodist(x=sdf2, measure = "haversine")
    sdf4 <- unique(as.vector(sdf3))
    sdf4 <- sdf4[sdf4 != 0]        # this is to remove the 0-distances 
    mean_dist <- mean(sdf4)
    distances <- c(distances, mean_dist)
    }  
    return(distances)
}

split(df1, df1$province) %>% give_distance()
#[1]  15379.21 793612.04

Например. функция даст вам средние значения расстояния для каждой провинции. Теперь мне не удалось заставить give_distance работать с sapply, но это уже должно быть более эффективно.

В отношении этого нить векторизованное решение вашей проблемы будет выглядеть следующим образом:

toCheck <- sapply(split(df1, df1$province), function(x){
                                            combn(rownames(x), 2, simplify = FALSE)})

names(toCheck) <- sapply(toCheck, paste, collapse = " - ")


sapply(toCheck, function(x){
               distm(df1[x[1],c("lon","lat")], df1[x[2],c("lon","lat")], 
                     fun = distHaversine)
                           })


  #    1 - 2      1 - 3      2 - 3      4 - 5      4 - 6      5 - 6 
  # 11429.10   22415.04   12293.48  634549.20 1188925.65  557361.28 

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

для вашего фактического набора данных toCheck станет вложенным списком, поэтому вам нужно настроить функцию, как показано ниже; Я не сделал toCheck имена чистыми для этого решения. (df2 можно найти в конце ответа).

df2 <- df2[order(df2$province),] #sorting may even improve performance
names(toCheck) <- paste("province", unique(df2$province))

toCheck <- sapply(split(df2, df2$province), function(x){
                                            combn(rownames(x), 2, simplify = FALSE)})

sapply(toCheck, function(x){ sapply(x, function(y){
  distm(df2[y[1],c("lon","lat")], df2[y[2],c("lon","lat")], fun = distHaversine)
})})

# $`province 1`
# [1]   11429.10   22415.04 1001964.84   12293.48 1013117.36 1024209.46
# 
# $`province 2`
# [1]  634549.2 1188925.7  557361.3
# 
# $`province 3`
# [1] 590083.2
# 
# $`province 4`
# [1] 557361.28 547589.19  11163.92

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

df2 <- data.frame(province = c(1, 1, 1, 2, 2, 2, 1, 3, 3, 4,4,4),
                  house = c(1, 2, 3, 4, 5, 6, 7, 10, 9, 8, 11, 12),
                  lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7, -85.6, -76.4, -75.4, -80.9, -85.7, -85.6), 
                  lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2, 40.1, 39.3, 60.8, 53.3, 40.2, 40.1))

Решение:

lapply(split(df1, df1$province), function(df){
  df <- Expand.Grid(df[, c("lat", "lon")], df[, c("lat", "lon")])
  mean(distHaversine(df[, 1:2], df[, 3:4]))
})

где Expand.Grid() взято из https://stackoverflow.com/a/30085602/3502164.

Объяснение:

1. Производительность

Я бы не стал использовать distm(), так как он преобразует функцию векторизованныйdistHaversine() в невекторизованную distm(). Если вы посмотрите на исходный код, вы увидите:

function (x, y, fun = distHaversine) 
{
   [...]
   for (i in 1:n) {
        dm[i, ] = fun(x[i, ], y)
    }
    return(dm)
}

В то время как distHaversine() отправляет «целый объект» в C, distm() отправляет данные «по строкам» в distHaversine() и, следовательно, заставляет distHaversine() делать то же самое при выполнении кода в C. Поэтому distm() не следует использовать. С точки зрения производительности я вижу больше вреда от использования функции-оболочки distm(), чем преимуществ.

2. Объяснение кода в «решении»:

а) Разделение на группы:

Вы хотите проанализировать данные по группе: провинция. Разделить на группы можно: split(df1, df1$province).

б) Группировка «сгустков столбцов»

Вы хотите найти все уникальные комбинации широты и долготы. Первое предположение может быть expand.grid(), но это не работает для нескольких столбцов. К счастью, мистер Флик позаботился об этом Функция expand.grid для data.frames в R.

Тогда у вас есть data.frame() всех возможных комбинаций, и вам просто нужно использовать mean(distHaversine(...)).

Моя первоначальная идея состояла в том, чтобы просмотреть исходный код distHaversine и воспроизвести его в функции, которую я буду использовать с proxy. Это будет работать следующим образом (обратите внимание, что lon должен быть первым столбцом):

library(geosphere)
library(dplyr)
library(proxy)

df1 <- data.frame(province = as.integer(c(1, 1, 1, 2, 2, 2)),
                  house = as.integer(c(1, 2, 3, 4, 5, 6)),
                  lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7), 
                  lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2))

custom_haversine <- function(x, y) {
  toRad <- pi / 180

  diff <- (y - x) * toRad
  dLon <- diff[1L]
  dLat <- diff[2L]

  a <- sin(dLat / 2) ^ 2 + cos(x[2L] * toRad) * cos(y[2L] * toRad) * sin(dLon / 2) ^ 2
  a <- min(a, 1)
  # return
  2 * atan2(sqrt(a), sqrt(1 - a)) * 6378137
}

pr_DB$set_entry(FUN=custom_haversine, names = "haversine", loop=TRUE, distance=TRUE)

average_dist <- df1 %>%
  select(-house) %>%
  group_by(province) %>%
  group_map(~ data.frame(avg=mean(proxy::dist(.x[ , c("lon", "lat")], method = "haversine"))))

Однако, если вы ожидаете миллионы строк на провинцию, proxy вероятно не сможет выделить промежуточные (нижние треугольные) матрицы. Поэтому я перенес код на C++ и в качестве бонуса добавил многопоточность:

РЕДАКТИРОВАТЬ: оказывается, помощник s2d был далек от оптимального, в этой версии теперь используются формулы, заданные здесь.

РЕДАКТИРОВАТЬ2: Я только что узнал о RcppThread, и его можно использовать для обнаружения пользовательского прерывания.

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel,RcppThread)]]

#include <cstddef> // size_t
#include <math.h> // sin, cos, sqrt, atan2, pow
#include <vector>

#include <RcppThread.h>
#include <Rcpp.h>
#include <RcppParallel.h>

using namespace std;
using namespace Rcpp;
using namespace RcppParallel;

// single to double indices for lower triangular of matrices without diagonal
void s2d(const size_t id, const size_t nrow, size_t& i, size_t& j) {
  j = nrow - 2 - static_cast<size_t>(sqrt(-8 * id + 4 * nrow * (nrow - 1) - 7) / 2 - 0.5);
  i = id + j + 1 - nrow * (nrow - 1) / 2 + (nrow - j) * ((nrow - j) - 1) / 2;
}

class HaversineCalculator : public Worker
{
public:
  HaversineCalculator(const NumericVector& lon,
                      const NumericVector& lat,
                      double& avg,
                      const int n)
    : lon_(lon)
    , lat_(lat)
    , avg_(avg)
    , n_(n)
    , cos_lat_(lon.length())
  {
    // terms for distance calculation
    for (size_t i = 0; i < cos_lat_.size(); i++) {
      cos_lat_[i] = cos(lat_[i] * 3.1415926535897 / 180);
    }
  }

  void operator()(size_t begin, size_t end) {
    // for Kahan summation
    double sum = 0;
    double c = 0;

    double to_rad = 3.1415926535897 / 180;

    size_t i, j;
    for (size_t ind = begin; ind < end; ind++) {
      if (RcppThread::isInterrupted(ind % static_cast<int>(1e5) == 0)) return;

      s2d(ind, lon_.length(), i, j);

      // haversine distance
      double d_lon = (lon_[j] - lon_[i]) * to_rad;
      double d_lat = (lat_[j] - lat_[i]) * to_rad;
      double d_hav = pow(sin(d_lat / 2), 2) + cos_lat_[i] * cos_lat_[j] * pow(sin(d_lon / 2), 2);
      if (d_hav > 1) d_hav = 1;
      d_hav = 2 * atan2(sqrt(d_hav), sqrt(1 - d_hav)) * 6378137;

      // the average part
      d_hav /= n_;

      // Kahan sum step
      double y = d_hav - c;
      double t = sum + y;
      c = (t - sum) - y;
      sum = t;
    }

    mutex_.lock();
    avg_ += sum;
    mutex_.unlock();
  }

private:
  const RVector<double> lon_;
  const RVector<double> lat_;
  double& avg_;
  const int n_;
  tthread::mutex mutex_;
  vector<double> cos_lat_;
};

// [[Rcpp::export]]
double avg_haversine(const DataFrame& input, const int nthreads) {
  NumericVector lon = input["lon"];
  NumericVector lat = input["lat"];

  double avg = 0;
  int size = lon.length() * (lon.length() - 1) / 2;
  HaversineCalculator hc(lon, lat, avg, size);

  int grain = size / nthreads / 10;
  RcppParallel::parallelFor(0, size, hc, grain);
  RcppThread::checkUserInterrupt();

  return avg;
}

Этот код не будет выделять промежуточную матрицу, он просто рассчитает расстояние для каждой пары того, что будет нижним треугольником, и в конце аккумулирует значения для среднего значения. См. здесь для части суммирования Кахана.

Если вы сохраните этот код, скажем, в haversine.cpp, то вы можете сделать следующее:

library(dplyr)
library(Rcpp)
library(RcppParallel)
library(RcppThread)

sourceCpp("haversine.cpp")

df1 %>%
  group_by(province) %>%
  group_map(~ data.frame(avg=avg_haversine(.x, parallel::detectCores())))
# A tibble: 2 x 2
# Groups:   province [2]
  province     avg
     <int>   <dbl>
1        1  15379.
2        2 793612.

Вот еще проверка на вменяемость:

pr_DB$set_entry(FUN=geosphere::distHaversine, names = "distHaversine", loop=TRUE, distance=TRUE)

df1 %>%
  select(-house) %>%
  group_by(province) %>%
  group_map(~ data.frame(avg=mean(proxy::dist(.x[ , c("lon", "lat")], method = "distHaversine"))))

Однако предостережение:

df <- data.frame(lon=runif (1e3, -90, 90), lat=runif (1e3, -90, 90))

system.time(proxy::dist(df, method = "distHaversine"))
   user  system elapsed 
 34.353   0.005  34.394

system.time(proxy::dist(df, method = "haversine"))
   user  system elapsed 
  0.789   0.020   0.809

system.time(avg_haversine(df, 4L))
   user  system elapsed 
  0.054   0.000   0.014

df <- data.frame(lon=runif (1e5, -90, 90), lat=runif (1e5, -90, 90))

system.time(avg_haversine(df, 4L))
   user  system elapsed 
 73.861   0.238  19.670

Вам, вероятно, придется подождать довольно долго, если у вас есть миллионы строк...

Я также должен упомянуть, что невозможно обнаружить прерывание пользователя внутри потоков, созданных через RcppParallel, поэтому, если вы начинаете расчет, вы должны либо дождаться его завершения, либо или полностью перезапустите R/RStudio. См. РЕДАКТИРОВАТЬ2 выше.

Относительно сложности

В зависимости от ваших фактических данных и количества ядер вашего компьютера, вы вполне можете в конечном итоге ждать несколько дней, пока расчет не будет завершен. Эта задача имеет квадратичную сложность (по провинции, так сказать). Эта строка:

int size = lon.length() * (lon.length() - 1) / 2;

означает количество (гаверсинус) вычислений расстояния, которые должны быть выполнены. Таким образом, если количество строк увеличивается в n раз, количество вычислений увеличивается в n^2 / 2 раз, грубо говоря.

Это невозможно оптимизировать; вы не можете рассчитать среднее N чисел без фактического вычисления каждого числа, и вам будет трудно найти что-то более быстрое, чем многопоточный код C++, так что вам придется либо подождать, или кинуть больше ядер на проблему, либо с одной машиной, либо со многими машинами, работающими вместе. Иначе эту проблему не решить.

Учитывая, что ваши данные содержат миллионы строк, это звучит как проблема «XY». т.е. ответ, который вам действительно нужен, не является ответом на вопрос, который вы задали.

Приведу аналогию: если вы хотите узнать среднюю высоту деревьев в лесу, вы не будете измерять каждое дерево. Вы просто измеряете достаточно большую выборку, чтобы убедиться, что ваша оценка имеет достаточно высокую вероятность быть настолько близкой к истинному среднему, насколько вам нужно.

Выполнение расчета методом грубой силы с использованием расстояния от каждого дома до каждого второго дома не только потребует чрезмерных ресурсов (даже с оптимизированным кодом), но и даст гораздо больше знаков после запятой, чем вам может понадобиться, или оправдано точностью данных (Координаты GPS обычно точны в лучшем случае с точностью до нескольких метров).

Поэтому я бы порекомендовал выполнять расчеты с размером выборки, который настолько велик, насколько это необходимо для уровня точности, которого требует ваша задача. Например, в следующем примере будет получена оценка для двух миллионов строк с точностью до 4 значащих цифр всего за несколько секунд. Вы можете повысить точность, увеличив размер выборки, но, учитывая неопределенность самих GPS-координат, я сомневаюсь, что это оправдано.

sample.size=1e6    
lapply(split(df1[3:4], df1$province), 
  function(x) {
    s1 = x[sample(nrow(x), sample.size, T), ]
    s2 = x[sample(nrow(x), sample.size, T), ]
    mean(distHaversine(s1, s2))
  })

Некоторые большие данные для тестирования:

N=1e6
df1 <- data.frame(
  province = c(rep(1,N),rep(2,N)),
  house = 1:(2*N),
  lat = c(rnorm(N,-76), rnorm(N,-85)), 
  lon = c(rnorm(N,39), rnorm(N,-55,2)))

Чтобы получить представление о точности этого метода, мы можем использовать начальную загрузку. Для следующей демонстрации я использую всего 100 000 строк данных, чтобы мы могли выполнить 1000 итераций начальной загрузки за короткое время:

N=1e5
df1 <- data.frame(lat = rnorm(N,-76,0.1), lon = rnorm(N,39,0.1))

dist.f = function(i) {
    s1 = df1[sample(N, replace = T), ]
    s2 = df1[sample(N, replace = T), ]
    mean(distHaversine(s1, s2))
    }

boot.dist = sapply(1:1000, dist.f)
mean(boot.dist)
# [1] 17580.63
sd(boot.dist)
# [1] 29.39302

hist(boot.dist, 20) 

т.е. для этих тестовых данных среднее расстояние составляет 17 580 +/- 29 м. Это коэффициент вариации 0,1%, что, вероятно, достаточно точно для большинства целей. Как я уже сказал, вы можете добиться большей точности, увеличив размер выборки, если вам это действительно нужно.

Вы можете использовать векторизованную версию расстояния гаверсинуса, например:

dist_haversine_for_dfs <- function (df_x, df_y, lat, r = 6378137) 
{
  if (!all(c("lat", "lon") %in% names(df_x))) {
    stop("parameter df_x does not have column 'lat' and 'lon'")
  }
  if (!all(c("lat", "lon") %in% names(df_y))) {
    stop("parameter df_x does not have column 'lat' and 'lon'")
  }
  toRad <- pi/180
  df_x <- df_x * toRad
  df_y <- df_y * toRad
  dLat <- df_y[["lat"]] - df_x[["lat"]]
  dLon <- df_y[["lon"]] - df_x[["lon"]]
  a <- sin(dLat/2) * sin(dLat/2) + cos(df_x[["lat"]]) * cos(df_y[["lat"]]) * 
    sin(dLon/2) * sin(dLon/2)
  a <- pmin(a, 1)
  dist <- 2 * atan2(sqrt(a), sqrt(1 - a)) * r
  return(dist)
}

Затем с помощью data.table и пакета arrangements (для более быстрой генерации комбинаций) можно сделать следующее:

library(data.table)
dt <- data.table(df1)
ids <- dt[, {
  comb_mat <- arrangements::combinations(x = house, k = 2)
  list(house_x = comb_mat[, 1],
       house_y = comb_mat[, 2])}, by = province]

jdt <- cbind(ids, 
             dt[ids$house_x, .(lon_x=lon, lat_x=lat)], 
             dt[ids$house_y, .(lon_y=lon, lat_y=lat)])

jdt[, dist := dist_haversine_for_dfs(df_x = jdt[, .(lon = lon.x, lat = lat.x)],
                                     df_y = jdt[, .(lon = lon.y, lat = lat.y)])]

jdt[, .(mean_dist = mean(dist)), by = province]

который выводит

   province mean_dist
1:        1  15379.21
2:        2 793612.04

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

library(data.table)
library(tidyverse)
library(spatialrisk)
library(optiRum)

# Expand grid
grid <- function(x){
  df <- x[, lat, lon]
  optiRum::CJ.dt(df, df)
}

Поскольку каждый элемент вывода представляет собой фрейм данных, purrr::map_dfr используется для связывания их вместе:

data.table(df1) %>%
  split(.$province) %>%
  map_dfr(grid, .id = "province") %>%
  mutate(distm = spatialrisk::haversine(lat, lon, i.lat, i.lon)) %>%
  filter(distm > 0) %>%
  group_by(province) %>%
  summarize(distm_mean = mean(distm))

Выход:

  province distm_mean
  <chr>         <dbl>
1 1            15379.
2 2           793612.