2016-06-10 156 views
1

有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) 

} 

回答

3

您可以使用sp-package中的spDists高效地完成此操作。是什麼讓它非常快速的是,實際的距離計算是在C中完成的,而不是在R中。解釋型語言(如R)中的循環非常緩慢。

1

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 
1

下面是一個使用data.tabledistGeo{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")