First, let's download some data similar to yours (I assume). This csv available online has almost 7,000 airports:
url <- "https://commondatastorage.googleapis.com/ckannet-storage/2012-07-09T214020/global_airports.csv"
library(RCurl)
txt <- getURL(url)
data <- read.csv(textConnection(txt), stringsAsFactors = FALSE)
For illustration purposes, we will use a small sample: the six airports in Jamaica.
dat <- subset(data, country == "Jamaica",
c("city", "country", "name", "latitude", "longitude"))
dat
# city country name latitude longitude
# 1745 Ocho Rios Jamaica Boscobel 18.40425 -76.96902
# 1746 Kingston Jamaica Norman Manley Intl 17.93567 -76.78750
# 1747 Montego Bay Jamaica Sangster Intl 18.50372 -77.91336
# 1748 Port Antonio Jamaica Ken Jones 18.19881 -76.53453
# 1749 Kingston Jamaica Tinson Pen 17.98856 -76.82376
# 5878 Negril Jamaica Negril Aerodrome 18.34000 -78.33556
Now let's have a look at your code. I will not review the math in earth.dist
, I'll assume it is correct. One beautiful thing about that function is that it is vectorized, i.e., you could give it n
-long vectors as inputs and it will compute n
distances in a single call. Unfortunately, the rest of your code does not take advantage of it. Instead, your double loop only calls earth.dist
with scalars at each time...
Instead of a double loop, you should be using the outer
function. Have a look at the doc (?outer
) if you are not familiar with it. The typical usage is outer(X, Y, FUN)
where X
and Y
are vectors and FUN
is a vectorized function. The output is a matrix Z
where Z[i, j]
is the result of FUN(X[i], Y[j])
. But what's brilliant about outer
is that it does not call FUN
as many times as there are entries in Z
(length(X) * length(Y)
). No, it calls is only once. How? Because FUN
is vectorized (a requirement) and outer
knows how to take advantage of it.
So, here is how we can massage your data a bit so we can use outer
. First, remember that outer
loops on the pairwise combinations from two vectors. In our case, we could use the names of the airports:
airport.names <- dat$name
so we will be calling outer(airport.names, airport.names, FUN = airport.dist)
. All that is left is to write airport.dist
: a vectorized function that will take as inputs two vectors of airport names and return their distances. We could first put the important data in a matrix with airport names as row names for easy access:
dat.mat <- as.matrix(dat[, c("latitude", "longitude")])
rownames(dat.mat) <- airport.names
Then define:
airport.dist <- function(name1, name2, data = dat.mat) {
lon1 <- data[name1, "longitude"]
lat1 <- data[name1, "latitude"]
lon2 <- data[name2, "longitude"]
lat2 <- data[name2, "latitude"]
return(earth.dist(lon1, lat1, lon2, lat2))
}
Then run outer
:
dist.mat <- outer(airport.names, airport.names, FUN = airport.dist)
and attach names to the columns and rows:
dimnames(dist.mat) <- list(airport.names, airport.names)
# Boscobel Norman Manley Intl Sangster Intl Ken Jones Tinson Pen Negril Aerodrome
# Boscobel 0.00000 55.58311 100.33081 51.30027 48.75737 144.54554
# Norman Manley Intl 55.58311 0.00000 134.79830 39.68379 7.02926 169.83797
# Sangster Intl 100.33081 134.79830 0.00000 149.58621 128.67955 48.17112
# Ken Jones 51.30027 39.68379 149.58621 0.00000 38.52863 191.03056
# Tinson Pen 48.75737 7.02926 128.67955 38.52863 0.00000 164.62138
# Negril Aerodrome 144.54554 169.83797 48.17112 191.03056 164.62138 0.00000
If you need to convince yourself that earth.dist
was called a single time, you could add a cat("HELLO\n")
somewhere inside its body (I did!). earth.dist
having been called only once, there is no need to say how fast the computation will be.
Finally, if you want to store the distances in a three column (airport1, airport2, distance) data.frame rather than a matrix, you can do:
d <- dist.mat
dist.df <- data.frame(airport1 = rownames(d)[row(d)],
airport2 = colnames(d)[col(d)],
distance = c(dist.mat))
I hope it helps! Don't hesitate to comment below if you have questions.