Hi Curt, 

Thanks for the help. 

According to that blogpost you sent, I am using the function gcd.hf using the Haversine formula. I wrapped it up in a function called CalcDists so that I can get a distance matrix between N sites. 

I don't know much about calculating distances, so does this seem like a good way to go (using Haversine that is)?


# 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
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
return(d) # Distance in km
}

# Fxn to calculate matrix of distances between each two sites
CalcDists <- function(latlongs) {
name <- list(rownames(latlongs), rownames(latlongs))
n <- nrow(latlongs)
z <- matrix(0, n, n, dimnames = name)
for (i in 1:n) {
for (j in 1:n) z[i, j] <- gcd.hf(long1 = latlongs[i, 1], 
lat1 = latlongs[i, 2], long2 = latlongs[j, 1], lat2 = latlongs[j,2])
}
z <- as.dist(z)
return(z)
}



Scott
On Monday, April 11, 2011 at 5:00 PM, Seeliger.Curt@epamail.epa.gov wrote: 
> > A comparison of some geographic distance calculations is provided at 
> > http://pineda-krch.com/2010/11/23/great-circle-distance-calculations-in-r/ 
> > , along with code for calculating the Vincenty inverse formula, which 
> > relies on the WGS-84 ellipsoid approximations.
> 
> You know, Scott, I should have included some test results of that method. Comparing the distances with Arc 9 indicates that the accuracy varies with location and whether there is a longitudinal difference in the two points. Comparing calculation results for points shifted 0 secs to 10 degrees North, West and Northwest from a 'base' point, the relative errors (defined as (Arc9.distance - Vincenty.distance)/Arc9.distance) range up to 0.08 in AK, AZ, CA, MT, NE, NM, UT, WA and WY, and range only up to 0.009 otherwise. In the special case of zero longitudinal offset (North-South distances only), the relative error ranges to 0.006 in those states and to 2E-7 otherwise. 
> 
> Let us know if you can do better, 
> cur
> 
>  -- 
>  Curt Seeliger, Data Ranger
>  Raytheon Information Services - Contractor to ORD
> seeliger.curt@epa.gov
>  541/754-4638
> 

	[[alternative HTML version deleted]]

