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

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

> Dat
  IndIDII      IndYear  WintLat  WintLong
1 BHS_265 BHS_265-2015 47.61025 -112.7210
2 BHS_265 BHS_265-2016 47.59884 -112.7089
3 BHS_770 BHS_770-2016 42.97379 -109.0400
4 BHS_770 BHS_770-2017 42.97129 -109.0367
5 BHS_770 BHS_770-2018 42.97244 -109.0509
6 BHS_377 BHS_377-2015 43.34744 -109.4821
7 BHS_377 BHS_377-2016 43.35559 -109.4445
8 BHS_377 BHS_377-2017 43.35195 -109.4566
9 BHS_377 BHS_377-2018 43.34765 -109.4892

Я хотел бы рассчитать евклидово расстояние между последовательными точками для каждого человека. Моя первоначальная задача заключалась в том, чтобы работать с dplyr, используя lead(), как показано ниже. Для функции distm требуется матрица, которую я не смог создать в dplyr. Можно ли сгенерировать и использовать матрицу в качестве аргумента для distm?

Dat %>% 
  group_by(IndIDII) %>% 
  mutate(WitnGeoDist = distm(as.matrix(c("WintLong", "WintLat")), lead(as.matrix(c("WintLong", "WintLat"))), fun = distVincentyEllipsoid))

В качестве альтернативы, есть ли другие возможности ...? Спасибо заранее.

Данные:

Dat <- structure(list(IndIDII = c("BHS_265", "BHS_265", "BHS_770", "BHS_770", 
"BHS_770", "BHS_377", "BHS_377", "BHS_377", "BHS_377"), IndYear = c("BHS_265-2015", 
"BHS_265-2016", "BHS_770-2016", "BHS_770-2017", "BHS_770-2018", 
"BHS_377-2015", "BHS_377-2016", "BHS_377-2017", "BHS_377-2018"
), WintLat = c(47.6102519805014, 47.5988417247191, 42.9737859090909, 
42.9712914772727, 42.9724390816327, 43.3474354347826, 43.3555934579439, 
43.3519543396226, 43.3476466990291), WintLong = c(-112.720994832869, 
-112.708887595506, -109.039964727273, -109.036693522727, -109.050923061224, 
-109.482114456522, -109.444522149533, -109.45659254717, -109.489241553398
)), class = "data.frame", row.names = c(NA, -9L))

r sp distance dplyr geosphere
person B. Davis    schedule 19.09.2018    source источник
comment
Если вы просто используете евклидово расстояние, вы можете использовать вычисление напрямую sqrt( (y2 - y1)^2 + (x2 - x1)^2 ), distm не требуется.   -  person SymbolixAU    schedule 20.09.2018

Ответы (2)


arrow_upward
1
arrow_downward

Вот другой метод, который group_by лучше использует и заставляет geosphere::distm работать, используя purrr::possibly. Это позволяет нам заполнить NA для строк, где расстояние не имеет смысла, потому что нет предыдущих значений для работы.

Dat <- structure(list(IndIDII = c("BHS_265", "BHS_265", "BHS_770", "BHS_770", "BHS_770", "BHS_377", "BHS_377", "BHS_377", "BHS_377"), IndYear = c("BHS_265-2015", "BHS_265-2016", "BHS_770-2016", "BHS_770-2017", "BHS_770-2018", "BHS_377-2015", "BHS_377-2016", "BHS_377-2017", "BHS_377-2018"), WintLat = c(47.6102519805014, 47.5988417247191, 42.9737859090909, 42.9712914772727, 42.9724390816327, 43.3474354347826, 43.3555934579439, 43.3519543396226, 43.3476466990291), WintLong = c(-112.720994832869, -112.708887595506, -109.039964727273, -109.036693522727, -109.050923061224, -109.482114456522, -109.444522149533, -109.45659254717, -109.489241553398)), class = "data.frame", row.names = c(NA, -9L))
library(tidyverse)
poss_dist <- possibly(geosphere::distm, otherwise = NA)
Dat %>%
  nest(WintLong, WintLat, .key = "coords") %>%
  group_by(IndIDII) %>%
  mutate(prev_coords = lag(coords)) %>%
  ungroup() %>%
  mutate(WitnGeoDist = map2_dbl(coords, prev_coords, poss_dist))
#> # A tibble: 9 x 5
#>   IndIDII IndYear      coords              prev_coords         WitnGeoDist
#>   <chr>   <chr>        <list>              <list>                    <dbl>
#> 1 BHS_265 BHS_265-2015 <data.frame [1 x 2~ <lgl [1]>                   NA 
#> 2 BHS_265 BHS_265-2016 <data.frame [1 x 2~ <data.frame [1 x 2~       1561.
#> 3 BHS_770 BHS_770-2016 <data.frame [1 x 2~ <lgl [1]>                   NA 
#> 4 BHS_770 BHS_770-2017 <data.frame [1 x 2~ <data.frame [1 x 2~        385.
#> 5 BHS_770 BHS_770-2018 <data.frame [1 x 2~ <data.frame [1 x 2~       1168.
#> 6 BHS_377 BHS_377-2015 <data.frame [1 x 2~ <lgl [1]>                   NA 
#> 7 BHS_377 BHS_377-2016 <data.frame [1 x 2~ <data.frame [1 x 2~       3180.
#> 8 BHS_377 BHS_377-2017 <data.frame [1 x 2~ <data.frame [1 x 2~       1059.
#> 9 BHS_377 BHS_377-2018 <data.frame [1 x 2~ <data.frame [1 x 2~       2690.

Создано 19 сентября 2018 г. пакетом REPEX (v0.2.0).

person Calum You    schedule 19.09.2018
comment
Я думаю, вы можете вызвать geosphere::distVincentyEllipsoid напрямую, а не проходить через distm (который создает матрицу расстояний). - person SymbolixAU; 20.09.2018

arrow_upward
1
arrow_downward

Вот подход sf и tidyverse, хотя я не думаю, что он самый чистый. Мне не удалось заставить geosphere::distm изящно обрабатывать отсутствующие значения (что позволило бы нам использовать group_by), поэтому вместо этого я прибег к использованию split с st_distance.

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

Dat <- structure(list(IndIDII = c("BHS_265", "BHS_265", "BHS_770", "BHS_770", "BHS_770", "BHS_377", "BHS_377", "BHS_377", "BHS_377"), IndYear = c("BHS_265-2015", "BHS_265-2016", "BHS_770-2016", "BHS_770-2017", "BHS_770-2018", "BHS_377-2015", "BHS_377-2016", "BHS_377-2017", "BHS_377-2018"), WintLat = c(47.6102519805014, 47.5988417247191, 42.9737859090909, 42.9712914772727, 42.9724390816327, 43.3474354347826, 43.3555934579439, 43.3519543396226, 43.3476466990291), WintLong = c(-112.720994832869, -112.708887595506, -109.039964727273, -109.036693522727, -109.050923061224, -109.482114456522, -109.444522149533, -109.45659254717, -109.489241553398)), class = "data.frame", row.names = c(NA, -9L))
library(tidyverse)
library(sf)

Dat %>%
  st_as_sf(coords = c("WintLong", "WintLat"), crs = 4326, remove = FALSE) %>%
  split(.$IndIDII) %>%
  map(function(df){
    dist <- st_distance(df[2:nrow(df), ], df[1:(nrow(df)- 1), ], by_element = TRUE)
    df %>% mutate(WitnGeoDist = c(NA, dist))
  }) %>%
  invoke(rbind, .x = .)
#> Simple feature collection with 9 features and 5 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -112.721 ymin: 42.97129 xmax: -109.0367 ymax: 47.61025
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#>           IndIDII      IndYear  WintLat  WintLong WitnGeoDist
#> BHS_265.1 BHS_265 BHS_265-2015 47.61025 -112.7210          NA
#> BHS_265.2 BHS_265 BHS_265-2016 47.59884 -112.7089   1561.4776
#> BHS_377.1 BHS_377 BHS_377-2015 43.34744 -109.4821          NA
#> BHS_377.2 BHS_377 BHS_377-2016 43.35559 -109.4445   3179.6929
#> BHS_377.3 BHS_377 BHS_377-2017 43.35195 -109.4566   1058.7986
#> BHS_377.4 BHS_377 BHS_377-2018 43.34765 -109.4892   2689.9938
#> BHS_770.1 BHS_770 BHS_770-2016 42.97379 -109.0400          NA
#> BHS_770.2 BHS_770 BHS_770-2017 42.97129 -109.0367    384.7117
#> BHS_770.3 BHS_770 BHS_770-2018 42.97244 -109.0509   1167.7996
#>                             geometry
#> BHS_265.1  POINT (-112.721 47.61025)
#> BHS_265.2 POINT (-112.7089 47.59884)
#> BHS_377.1 POINT (-109.4821 43.34744)
#> BHS_377.2 POINT (-109.4445 43.35559)
#> BHS_377.3 POINT (-109.4566 43.35195)
#> BHS_377.4 POINT (-109.4892 43.34765)
#> BHS_770.1   POINT (-109.04 42.97379)
#> BHS_770.2 POINT (-109.0367 42.97129)
#> BHS_770.3 POINT (-109.0509 42.97244)

Создано 19 сентября 2018 г. пакетом REPEX (v0.2.0).

person Calum You    schedule 19.09.2018