分组数据中点间欧氏距离的计算
分组数据中点间欧氏距离的计算
提问于 2018-09-19 22:14:28
在下面的数据中(包括在dput中)，我对三个人(IndIDII)重复观察(lat和long)。注意，每个人都有不同的位置。
> 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))回答 2
Stack Overflow用户
回答已采纳
发布于 2018-09-19 23:06:05
下面是一种不同的方法，它更好地利用了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.由reprex封装创建于2018-09-19 (v0.2.0)。
Stack Overflow用户
发布于 2018-09-19 22:44:20
这是一种sf和tidyverse方法，尽管我不认为它是最干净的。我无法让geosphere::distm优雅地处理丢失的值(这会让我们使用group_by)，所以我转而使用split和st_distance。
这些步骤基本上是将坐标转换为点几何图形，在分组列上拆分以生成数据列表，使用添加距离列的函数在此列表中映射，然后将数据重新组合在一起。
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)由reprex封装创建于2018-09-19 (v0.2.0)。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/52414879
复制相关文章
相似问题
腾讯云开发者
