在广泛的数据框中计算每对坐标之间的距离。

Question

在广泛的数据框中计算每对坐标之间的距离。

5

我想计算两个空间坐标集合（在我的虚假数据集中是program和admin）之间的距离。数据以宽格式存储，因此两对坐标在同一行。

library(sp)
set.seed(1)
n <- 100
program.id <- seq(1, n)
c1 <- cbind(runif(n, -90, 90), runif(n, -180, 180))
c2 <- cbind(runif(n, -90, 90), runif(n, -180, 180))
dat <- data.frame(cbind(program.id, c1, c2))
names(dat) <- c("program.id", "program.lat", "program.long", "admin.lat", "admin.long")
head(dat)
#       program.id program.lat program.long  admin.lat  admin.long
# 1              1   -42.20844     55.70061 -41.848523   62.536404
# 2              2   -23.01770    -52.84898 -50.643849 -145.851172
# 3              3    13.11361    -82.70635   3.023431   -2.665397
# 4              4    73.47740    177.36626 -41.588893  -13.841337
# 5              5   -53.69725     48.05758 -57.389701  -44.922049
# 6              6    71.71014   -103.24507   3.343705  176.795719

我知道如何使用 sp 包创建 program 或 admin 之间的距离矩阵：

ll <- c("program.lat", "program.long")
coords <- dat[ll]
dist <- apply(coords, 1, 
              function(eachPoint) spDistsN1(as.matrix(coords),
                                            eachPoint, longlat=TRUE))

但我想要做的是创建一个nx1的距离向量(dist.km)，其中包含每对坐标之间的距离，并将其添加到dat中。

#       program.id program.lat program.long  admin.lat  admin.long  dist.km
# 1              1   -42.20844     55.70061 -41.848523   62.536404   567.35
# 2              2   -23.01770    -52.84898 -50.643849 -145.851172  8267.86
# ...

有什么建议吗？我花了一段时间查看了旧的SO问题，但似乎没有完全正确的答案。如果有更好的解决方法，请告诉我。

更新

@Amit的解决方案适用于我的玩具数据集：

apply(dat,1,function(x) spDistsN1(matrix(x[2:3],nrow=1),x[3:4],longlat=TRUE))

我认为我需要交换经纬度的顺序，使经度在纬度前面。这可以通过 ?spDistsN1 实现：

pts: A matrix of 2D points, first column x/longitude, second column y/latitude, or a SpatialPoints or SpatialPointsDataFrame object

此外，除非我误解了逻辑，否则我认为Amit的解决方案应该获取第2到3列和第4到5列，而不是第2到3列和第3到4列。

我的挑战现在是将此应用于我的实际数据。下面是部分复制的内容。

library(sp)
dat <- structure(list(ID = 1:4, 
                      subcounty = c("a", "b", "c", "d"), 
                      pro.long = c(33.47627919, 31.73605491, 31.54073482, 31.51748984), 
                      pro.lat = c(2.73996953, 3.26530095, 3.21327597, 3.17784981), 
                      sub.long = c(33.47552, 31.78307, 31.53083, 31.53083), 
                      sub.lat = c(2.740362, 3.391209, 3.208736, 3.208736)), 
                 .Names = c("ID", "subcounty", "pro.long", "pro.lat", "sub.long", "sub.lat"),     
                 row.names = c(NA, 4L), class = "data.frame")
head(dat) 
#     ID subcounty pro.long  pro.lat sub.long  sub.lat
#  1   1         a 33.47628 2.739970 33.47552 2.740362
#  2   2         b 31.73605 3.265301 31.78307 3.391209
#  3   3         c 31.54073 3.213276 31.53083 3.208736
#  4   4         d 31.51749 3.177850 31.53083 3.208736
apply(dat, 1, function(x) spDistsN1(matrix(x[3:4], nrow=1),
                                    x[5:6],
                                    longlat=TRUE))

我遇到了这个错误：Error in spDistsN1(matrix(x[3:4], nrow = 1), x[5:6], longlat = TRUE) : pts必须是数字

我很困惑，因为这些列都是数字类型的：

> is.numeric(dat$pro.long)
[1] TRUE
> is.numeric(dat$pro.lat)
[1] TRUE
> is.numeric(dat$sub.long)
[1] TRUE
> is.numeric(dat$sub.lat)
[1] TRUE

- Eric Green

1

你尝试过这个吗：apply(dat,1,function(x) spDistsN1(matrix(x[2:3],nrow=1),x[3:4],longlat=TRUE))？ - amit

@amit，我没有。我想答案可能涉及到其中一个apply函数，但我不知道矩阵的正确规范。这似乎是解决方案。如果你想添加一个答案，我很乐意接受它。 - Eric Green

只要它能够正常工作并有帮助，我就很满意。我不太在意声誉这种东西，但还是谢谢你的提供。 - amit

我相信根据 sp 帮助文档，我需要更改我的经度和纬度列的顺序："pts 一个二维点矩阵，第一列为 x/经度，第二列为 y/纬度"。 - Eric Green

@amit，我想你的意思是要获取第2列到第3列和第4列到第5列，而不是第3列到第4列。对吗？ - Eric Green

2个回答

4

使用 data.table 和 geosphere 可以提供更快的解决方案。

library(data.table)
library(geosphere)

setDT(dat)[ , dist_km := distGeo(matrix(c(pro.long, pro.lat), ncol = 2), 
                                  matrix(c(sub.long, sub.lat), ncol = 2))/1000]

基准测试:

library(sp)

jlhoward <- function(dat) { dat$dist.km <- sapply(1:nrow(dat),function(i)
                             spDistsN1(as.matrix(dat[i,3:4]),as.matrix(dat[i,5:6]),longlat=T)) }

rafa.pereira <- function(dat2) { setDT(dat2)[ , dist_km := distGeo(matrix(c(pro.long, pro.lat), ncol = 2), 
                                                                 matrix(c(sub.long, sub.lat), ncol = 2))/1000] }


> system.time( jlhoward(dat) )
   user  system elapsed 
   8.94    0.00    8.94 

> system.time( rafa.pereira(dat) )
   user  system elapsed 
   0.07    0.00    0.08

数据

dat <- structure(list(ID = 1:4, 
                      subcounty = c("a", "b", "c", "d"), 
                      pro.long = c(33.47627919, 31.73605491, 31.54073482, 31.51748984), 
                      pro.lat = c(2.73996953, 3.26530095, 3.21327597, 3.17784981), 
                      sub.long = c(33.47552, 31.78307, 31.53083, 31.53083), 
                      sub.lat = c(2.740362, 3.391209, 3.208736, 3.208736)), 
                 .Names = c("ID", "subcounty", "pro.long", "pro.lat", "sub.long", "sub.lat"),     
                 row.names = c(NA, 4L), class = "data.frame")

# enlarge dataset to 40,000 pairs
dat <- dat[rep(seq_len(nrow(dat)), 10000), ]

- rafa.pereira

1

Rafa，谢谢你的信息和回答。你的解决方案肯定更快！ - jazzurro

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jlhoward · Accepted Answer

你遇到的问题是apply(...)将第一个参数强制转换为矩阵。按照定义，矩阵必须具有相同数据类型的所有元素。由于dat中的一列(dat$subcounty)是字符型，apply(...)会将所有东西都转换为字符型。在你的测试数据集中，所有数据都是数值型，所以你没有遇到这个问题。

以下代码应该有效：

dat$dist.km <- sapply(1:nrow(dat),function(i)
                spDistsN1(as.matrix(dat[i,3:4]),as.matrix(dat[i,5:6]),longlat=T))