使用简单的for循环处理空间数据

Question

很抱歉，这是一个for循环101的问题。我无法编写一个简单的for循环来生成基于经纬度数据的城市之间距离的表格。

locations <-read.csv("distances.csv")

locations返回以下表格：

       City Type       long      lat
1 Sheffield  EUR  -1.470085 53.38113
2        HK WRLD 114.109497 22.39643
3    Venice  EUR  12.315515 45.44085
4  New York WRLD -74.005941 40.71278

在此任务的特定部分中，我要制作一张城市之间距离（以千米为单位）的表格，并呈现成相关矩阵的形式，其中对角线上的数值都是0（即所有城市到自身的距离均为0）。

为此，我使用了sp包，该包需要一个经度-纬度数值的矩阵，所以我可以按以下方法删除文本：

datmax <- data.matrix(locations)
datmax2 <- datmax[,-1:-2]

spDistsN1工具可以通过比较矩阵中所有城市到一个特定城市的距离，来获取这些信息。显然，我可以使用下面的表达式来获取所有城市到谢菲尔德（城市或行号#1）的距离：

km <- spDistsN1(datmax2, datmax2[1,], longlat=TRUE)

这将正确给出：

[1]    0.000 9591.009 1329.882 5436.133

然而，为了达到我想要的相关矩阵样式输出，我希望对每个城市都进行此操作。因此，我尝试编写了一个for循环：

for (i in 1:nrow(datmax2)){
  kmnew <- spDistsN1(datmax2, datmax2[i,], longlat=TRUE)
}

这为我提供了纽约的正确数值：

[1]  5436.133 12967.023  6697.541     0.000

我猜测在循环中我把一个城市写成了另一个城市。感谢您帮助我找到错误所在，非常感谢。

- RichS

Answer 1

首先声明一个矩阵，并使用迭代器 i 来指示要填充的行：

kmnew <- matrix(NA, nrow=4, ncol=4)
for (i in 1:nrow(datmax2)){
  kmnew[i,] <- spDistsN1(datmax2, datmax2[i,], longlat=TRUE)
}

colnames(kmnew) <- locations$City
rownames(kmnew) <- locations$City

结果

> kmnew

          Sheffield        HK   Venice  New York
Sheffield     0.000  9591.009 1329.882  5436.134
HK         9591.009     0.000 9134.698 12967.024
Venice     1329.882  9134.698    0.000  6697.541
New York   5436.134 12967.024 6697.541     0.000

- Dominic Comtois

Answer 2

我不确定这是否是您正在寻找的内容

。

library(sp)

# Provide data for reproducibility
locations <- data.frame(City=c("Sheffield", "HK", "Venice", "New York"),
                    Type=c("EUR", "WRLD", "EUR", "WRLD"),
                    long=c(-1.470085, 114.109497, 12.315515, -74.005941),
                    lat=c(53.38113, 22.39643, 45.44085, 40.71278))

km <- apply(as.matrix(locations[, c(-1, -2)]), 1, function(x){
  spDistsN1(as.matrix(locations[, c(-1, -2)]), x, longlat=TRUE)
})

km <- data.frame(locations[, 1],  km)
names(km) <- c("City", as.character(locations[, 1]))
km

结果

       City Sheffield        HK   Venice  New York
1 Sheffield     0.000  9591.009 1329.882  5436.134
2        HK  9591.009     0.000 9134.698 12967.024
3    Venice  1329.882  9134.698    0.000  6697.541
4  New York  5436.134 12967.024 6697.541     0.000

- dimitris_ps

---

这里的所有答案都非常有用。Dominic对我来说是最容易掌握的，因为它基于我的现有思维过程。"apply"是一个优雅的解决方案，感谢Dimitris，而Nicola我会查看这个包，尽管你可能能够从我的尝试中看出来，超越二维思维可能需要我一些时间！ - RichS

Answer 3

你可以尝试使用geosphere包中的distm函数：

 distm(datmax2)
 #        [,1]     [,2]    [,3]     [,4]
 #[1,]       0  9586671 1329405  5427956
 #[2,] 9586671        0 9130036 12962132
 #[3,] 1329405  9130036       0  6687416
 #[4,] 5427956 12962132 6687416        0

它返回以米为单位的距离，并考虑地球的几何形状。

- nicola

---

是的，但它假设地球是一个球体；sp::spDistsN1则假设它是一个椭球体（WGS84）。 - Edzer Pebesma