测量坐标的纬度和经度对之间的距离[重复]答案

【问题标题】：Measure distance between pairs of latitude and longitude of co-ordinates [duplicate]测量坐标的纬度和经度对之间的距离[重复]
【发布时间】：2016-10-11 02:52:00
【问题描述】：

有 43000 对坐标，正如您可以想象的那样，使用下面的嵌套循环运行这是一个非常漫长的过程，我只是不知道如何使这项工作更快并排除循环。

以下是我的代码，它基于我在网上找到的内容以及我对使用 R 的基础知识非常古老且有限的知识。请教我如何在更高效、更优雅的代码块中完成此操作，以帮助我完全有可能。

for (i in 1:length(FixedAssets$Lat)) {
    for (p  in 1:length(FixedAssets$Lat)) {
        dist <- distCosine(as.vector(c(FixedAssets$Longitude[i],FixedAssets$Latitude[i]), mode = "any"), as.vector(c(FixedAssets$Longitude[p],FixedAssets$Latitude[p]), mode = "any"), r = 6378137)
    if (dist < 200) { FixedAssets$prox200[i] = FixedAssets$prox200[i] + 1 }
    if (dist < 500) { FixedAssets$prox500[i] = FixedAssets$prox500[i] + 1 }
    if (dist < 1000) { FixedAssets$prox1000[i] = FixedAssets$prox1000[i] + 1 }
    if (dist < 2000) { FixedAssets$prox2000[i] = FixedAssets$prox2000[i] + 1 }
    if (dist < 3000) { FixedAssets$prox3000[i] = FixedAssets$prox3000[i] + 1 }
    if (dist < 5000) { FixedAssets$prox5000[i] = FixedAssets$prox5000[i] + 1 }
    if (dist < 10000) { FixedAssets$prox10000[i] = FixedAssets$prox10000[i] + 1 }
    if (dist < 20000) { FixedAssets$prox20000[i] = FixedAssets$prox20000[i] + 1 }    
    }
print(i)

}

【问题讨论】：

标签： r gis distance

【解决方案1】：

您可以使用 sp-package 中的spDists 有效地执行此操作。它之所以非常快，是因为实际距离计算是在 C 中完成的，而不是在 R 中。R 等解释语言中的循环非常慢。

【讨论】：

【解决方案2】：

spDist 的替代品可以是来自geoshere 包的distm。但是，计算距离往往需要更多时间：

library(maps)
library(geosphere)
library(sp)

data(world.cities)
system.time(d <- distm(world.cities[1:10000, c("long", "lat")]))
#  User      System verstrichen 
# 69.72        0.59       70.77 

system.time(d2 <- spDists(as.matrix(world.cities[1:10000, c("long", "lat")]), longlat = T))
# User      System verstrichen 
# 66.11        0.48       66.89 

as.dist(round(d[1:5, 1:5]/1000, 1))
#        1      2      3      4
# 2    1.5                     
# 3 3589.8 3588.7              
# 4 1327.5 1326.6 2326.7       
# 5  105.9  104.4 3510.2 1270.1
as.dist(round(d2[1:5, 1:5],1))
#        1      2      3      4
# 2    1.5                     
# 3 3592.9 3591.8              
# 4 1328.4 1327.5 2328.6       
# 5  105.6  104.1 3513.3 1270.8

【讨论】：

【解决方案3】：

这是一个非常快速的解决方案，它使用data.table 和distGeo{geosphere} 来计算椭圆体上的距离。

library(geosphere)
library(data.table)

setDT(dt)[ , dist_km := distGeo(matrix(c(long, lat), ncol = 2), 
                                  matrix(c(long_dest, lat_dest), ncol = 2))/1000]

可重现示例的数据

library(maps)
library(reshape)  

# Load world cities data and keep only 300 cities, which will give us 90,000 pairs
data(world.cities)
world.cities <- world.cities[1:300,] 

# create all possible pairs of origin-destination in a long format
dt <- expand.grid.df(world.cities,world.cities)
names(dt)[10:11] <- c("lat_dest","long_dest")

【讨论】：