我想计算两个链接的空间坐标集之间的距离（假冒数据集中的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

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

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 
# ...

有什么建议吗？我花了一段时间去解决老问题，但似乎没有什么正确的。很高兴被证明是错误的。

更新

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

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

但我想我需要交换LAT的顺序，LAT长柱这么久的长期订单在拉特之前。从?spDistsN1：

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

而且，除非我误解的逻辑，我认为Amit的解决方案应该抓住的cols [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 must be numeric

我很困惑，因为这些列数字：

> 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

来源

2014-03-04 Eric Green

你尝试：应用（DAT，1，函数（x）的spDistsN1（矩阵（X [2：3]， nrow = 1），x [3：4]，longlat = TRUE））？ – amit

@amit，我没有。我认为答案可能涉及应用功能之一，但我不知道矩阵的正确规格。这看起来是解决方案。如果您想添加答案，我很乐意接受它。 –

只要它有用并且有帮助 - 我很高兴。我不太在意这种声誉，但感谢提供。 – amit

您遇到的问题是，apply(...)胁迫的第一个参数矩阵。根据定义，矩阵必须包含相同数据类型的所有元素。由于dat（dat$subcounty）中的一列是char，因此apply(...)将所有内容强制转换为char。在你的测试数据集中，一切都是数字的，所以你没有这个问题。

这应该工作：

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

来源

2014-03-05 08:35:32 jlhoward

感谢您解释@jlhoward。这工作。非常感激。 –

因为我有类似的情况，所以今天我遇到了这个解决方案。我喜欢这个想法。我想知道我们是否可以让它工作得更好。我有一个像2GB这样的大数据集，并且用data.table尝试了这段代码。处理过程实际上已经进行了一段时间。对于每一行，我们要求R创建两个矩阵并处理计算。我宁愿认为创建SPDF并处理相同的工作。至少对于每一行，我们不必将DF转换为矩阵。任何想法？我也想知道是否有另一个功能更快地处理同一个工作。 – jazzurro

@jazzurro，我相信有一个更快的解决方案，使用'data.table'和'geosphere' http://stackoverflow.com/questions/36817423/how-to-efficiently-calculate-distance-between-pair-of-coordinates -using-data-tab –

有使用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), ]

来源

2016-06-10 11:50:29

拉法，感谢您的留言和答案。你的解决方案肯定会更快！ – jazzurro

计算宽数据帧中每对坐标之间的距离

回答

数据

相关问题