get a matrix of pairwise geographic distances between locations

CalcDists.R

# Convert degrees to radians
deg2rad <- function(deg) return(deg*pi/180)

# Calculates the geodesic distance between two points specified by 
# radian latitude/longitude using the Haversine formula
# Ouputs distance between sites 1 and 2 as meters
gcd.hf <- function(long1, lat1, long2, lat2) {
  R <- 6371 # Earth mean radius [km]
  delta.long <- (long2 - long1)
  delta.lat <- (lat2 - lat1)
  a <- sin(delta.lat/2)^2 + cos(lat1) * cos(lat2) * sin(delta.long/2)^2
  c <- 2 * asin(min(1,sqrt(a)))
  d = (R * c)*1000
  return(d) # Distance in meters
}

# Fxn to calculate matrix of distances between each two sites
# INPUT: a data frame in which longs are in first column and lats in second column
# OUTPUT: a distance matrix (class dist) between all pairwise sites
# Output distances are in meters
CalcDists <- function(longlats) {
	name <- list(rownames(longlats), rownames(longlats))
	n <- nrow(longlats)
	z <- matrix(0, n, n, dimnames = name)
    for (i in 1:n) {
        for (j in 1:n) z[i, j] <- gcd.hf(long1 = deg2rad(longlats[i, 2]), 
            lat1 = deg2rad(longlats[i, 1]), long2 = deg2rad(longlats[j, 2]), 
            lat2 = deg2rad(longlats[j, 1]))
    }
	z <- as.dist(z)
    return(z)
}

# E.g.
longlats <- data.frame(long = rnorm(10), lat = rnorm(10))
dists <- CalcDists(longlats)

	# convert to kilometers from meters
dists/1000