geosphere/R:计算用4个点定义的两个大圆的交点
geosphere/R: Calculate intersection points between two great circles defined with 4 points
我正在尝试使用 geosphere
来计算以数据帧格式给出的两个大圆之间的交点,如下所示:
library(dplyr)
library(geosphere)
df <- data.frame(
# first line
ln1_lonb = 1:4,
ln1_lone = 2:5,
ln1_latb = 10:13,
ln1_late = 11:14,
# second line
ln2_lonb = 5:8,
ln2_lone = 6:9,
ln2_latb = 14:17,
ln2_late = 15:18
)
我尝试使用 geosphere
中的 gcintersect
函数,该函数将矩阵作为输入。为了在数据框中使用它,我使用了 cbind
,但似乎 mutate 不适用于此:
df %>%
mutate(
int_points = gcIntersect(
cbind(ln1_lonb, ln1_latb),
cbind(ln1_lone, ln1_late),
cbind(ln2_lonb, ln2_latb),
cbind(ln2_lone, ln2_late)
)
)
>Error: Column `int_points` must be length 4 (the number of rows) or one, not 16
错误可能是由于输出比预期的要长(数据帧的行数)。所以我试着把它放在一个列表中:
df %>%
mutate(
int_points = list(gcIntersect(
cbind(ln1_lonb, ln1_latb),
cbind(ln1_lone, ln1_late),
cbind(ln2_lonb, ln2_latb),
cbind(ln2_lone, ln2_late)
))
)
在这里我再次看到输出是所有组合,而不是获取每行 2 个交点的 4 个坐标。
预期输出将是新列中单元格中的列表,其中包含两个点的坐标。
是否有不使用循环(或purrr
)的解决方案,因为这会比mutate
慢得多。
第一行中 int_points 的值应与此输出相同:
gcIntersect(cbind(1,2), cbind(10,11), cbind(5,6), cbind(14,15))
我们可以做一个rowwise
library(tidyverse)
library(geosphere)
df %>%
rowwise %>%
mutate( int_points = list(gcIntersect(
cbind(ln1_lonb, ln1_latb),
cbind(ln1_lone, ln1_late),
cbind(ln2_lonb, ln2_latb),
cbind(ln2_lone, ln2_late)
))) %>%
ungroup %>%
mutate(int_points = map(int_points, as_tibble)) %>%
unnest(int_points)
# A tibble: 4 x 12
# ln1_lonb ln1_lone ln1_latb ln1_late ln2_lonb ln2_lone ln2_latb ln2_late lon1 lat1 lon2 lat2
# <int> <int> <int> <int> <int> <int> <int> <int> <dbl> <dbl> <dbl> <dbl>
#1 1 2 10 11 5 6 14 15 -176. -12.6 3.61 12.6
#2 2 3 11 12 6 7 15 16 -175. -13.6 4.60 13.6
#3 3 4 12 13 7 8 16 17 -174. -14.6 5.60 14.6
#4 4 5 13 14 8 9 17 18 -173. -15.6 6.59 15.6
考虑到地圈方法是矢量化的(正如人们在 R 中所期望的那样),我会这样做
gcIntersect(df[,c(1,3)], df[,c(2,4)], df[,c(5,7)], df[,c(6,8)])
# lon1 lat1 lon2 lat2
#[1,] -176.3902 -12.58846 3.609757 12.58846
#[2,] -175.3968 -13.58017 4.603188 13.58017
#[3,] -174.4023 -14.57291 5.597652 14.57291
#[4,] -173.4070 -15.56648 6.592953 15.56648
你也可以先这样重组
d <- df[,c("ln1_lonb", "ln1_latb", "ln1_lone", "ln1_late", "ln2_lonb", "ln2_latb", "ln2_lone", "ln2_late")]
gcIntersect(d[,c(1,2)], d[,c(3,4)], d[,c(5,6)], d[,c(7:8)])
或者像这样
begin1 <- df[,c("ln1_lonb", "ln1_latb")]
end1 <- df[,c("ln1_lone", "ln1_late")]
begin2 <- df[,c("ln2_lonb", "ln2_latb")]
end2 <- df[,c("ln2_lone", "ln2_late")]
gcIntersect(begin1, end1, begin2, end2)
我知道你特别要求 mutate
的解决方案,但我为其他想要清晰简单的解决方案的人添加了这些解决方案。
我正在尝试使用 geosphere
来计算以数据帧格式给出的两个大圆之间的交点,如下所示:
library(dplyr)
library(geosphere)
df <- data.frame(
# first line
ln1_lonb = 1:4,
ln1_lone = 2:5,
ln1_latb = 10:13,
ln1_late = 11:14,
# second line
ln2_lonb = 5:8,
ln2_lone = 6:9,
ln2_latb = 14:17,
ln2_late = 15:18
)
我尝试使用 geosphere
中的 gcintersect
函数,该函数将矩阵作为输入。为了在数据框中使用它,我使用了 cbind
,但似乎 mutate 不适用于此:
df %>%
mutate(
int_points = gcIntersect(
cbind(ln1_lonb, ln1_latb),
cbind(ln1_lone, ln1_late),
cbind(ln2_lonb, ln2_latb),
cbind(ln2_lone, ln2_late)
)
)
>Error: Column `int_points` must be length 4 (the number of rows) or one, not 16
错误可能是由于输出比预期的要长(数据帧的行数)。所以我试着把它放在一个列表中:
df %>%
mutate(
int_points = list(gcIntersect(
cbind(ln1_lonb, ln1_latb),
cbind(ln1_lone, ln1_late),
cbind(ln2_lonb, ln2_latb),
cbind(ln2_lone, ln2_late)
))
)
在这里我再次看到输出是所有组合,而不是获取每行 2 个交点的 4 个坐标。
预期输出将是新列中单元格中的列表,其中包含两个点的坐标。
是否有不使用循环(或purrr
)的解决方案,因为这会比mutate
慢得多。
第一行中 int_points 的值应与此输出相同:
gcIntersect(cbind(1,2), cbind(10,11), cbind(5,6), cbind(14,15))
我们可以做一个rowwise
library(tidyverse)
library(geosphere)
df %>%
rowwise %>%
mutate( int_points = list(gcIntersect(
cbind(ln1_lonb, ln1_latb),
cbind(ln1_lone, ln1_late),
cbind(ln2_lonb, ln2_latb),
cbind(ln2_lone, ln2_late)
))) %>%
ungroup %>%
mutate(int_points = map(int_points, as_tibble)) %>%
unnest(int_points)
# A tibble: 4 x 12
# ln1_lonb ln1_lone ln1_latb ln1_late ln2_lonb ln2_lone ln2_latb ln2_late lon1 lat1 lon2 lat2
# <int> <int> <int> <int> <int> <int> <int> <int> <dbl> <dbl> <dbl> <dbl>
#1 1 2 10 11 5 6 14 15 -176. -12.6 3.61 12.6
#2 2 3 11 12 6 7 15 16 -175. -13.6 4.60 13.6
#3 3 4 12 13 7 8 16 17 -174. -14.6 5.60 14.6
#4 4 5 13 14 8 9 17 18 -173. -15.6 6.59 15.6
考虑到地圈方法是矢量化的(正如人们在 R 中所期望的那样),我会这样做
gcIntersect(df[,c(1,3)], df[,c(2,4)], df[,c(5,7)], df[,c(6,8)])
# lon1 lat1 lon2 lat2
#[1,] -176.3902 -12.58846 3.609757 12.58846
#[2,] -175.3968 -13.58017 4.603188 13.58017
#[3,] -174.4023 -14.57291 5.597652 14.57291
#[4,] -173.4070 -15.56648 6.592953 15.56648
你也可以先这样重组
d <- df[,c("ln1_lonb", "ln1_latb", "ln1_lone", "ln1_late", "ln2_lonb", "ln2_latb", "ln2_lone", "ln2_late")]
gcIntersect(d[,c(1,2)], d[,c(3,4)], d[,c(5,6)], d[,c(7:8)])
或者像这样
begin1 <- df[,c("ln1_lonb", "ln1_latb")]
end1 <- df[,c("ln1_lone", "ln1_late")]
begin2 <- df[,c("ln2_lonb", "ln2_latb")]
end2 <- df[,c("ln2_lone", "ln2_late")]
gcIntersect(begin1, end1, begin2, end2)
我知道你特别要求 mutate
的解决方案,但我为其他想要清晰简单的解决方案的人添加了这些解决方案。