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.c...@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.c...@epa.gov
> 541/754-4638
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
