请看下面的小工作示例:
#### Pseudo data
nobs1 <- 4000
nobs2 <- 5000
mylon1 <- runif(nobs1, min=0, max=1)-76
mylat1 <- runif(nobs1, min=0, max=1)+37
mylon2 <- runif(nobs2, min=0, max=1)-76
mylat2 <- runif(nobs2, min=0, max=1)+37
#### define a distance function
thedistance <- function(lon1, lat1, lon2, lat2) {
R <- 6371 # Earth mean radius [km]
delta.lon <- (lon2 - lon1)
delta.lat <- (lat2 - lat1)
a <- sin(delta.lat/2)^2 + cos(lat1) * cos(lat2) * sin(delta.lon/2)^2
c <- 2 * asin(min(1,sqrt(a)))
d = R * c
return(d) …