Answer to great circle distance and find.neighbor

[Jason Pither <pitherj@biology.queensu.ca>]

Hi:

Thanks to Dr. M Sawada for providing a detailed response (below) to my
query (below), and others who provided various functions to calculate
the Great Circle Distance in S (below).

MY TAKE ON RESPONSES:

Because I will have to map my results at some point, I need to use an
appropriate projection that won't distort areas too much at a
continental scale (an alternative projection would need to be used for
accurate representation of distances).  The projection provided by Dr.
Sawada does this (below).

If, however, you are interested only in finding proper neighbor
distances using lat / longs that are spread over a large area on the
globe and not actually worrying about projecting results on a map, then
it would still be useful to find out how to replace the internal C
algorithms (findDist, findAllDist) in the "find.neighbor" function to
accurately calculate neighbor distances using the great circle distance.

Thanks,
Jason.

*****
ORIGINAL QUERY:

I'm using S-Plus Spatial Stats V1.5 for Windows NT.

I would like the "find.neighbor" function to calculate
great-circle-distances (using decimal degrees) rather than cartesian
distances.

I know that there are map projections that I could use to project my
decimal degrees, but I'm doing continental-scale analyses, and
distortions in distances will occur at this scale with any single
projection (I think).

I know that the find.neighbor function calls on a C function
"findAllDist".  This is presumably where the alteration needs to occur.
***********

Dr. Sawada's Great Circle Distance function:

great.circle.distance.f<-function(cLonLat, oLonLat)
{
#cLonLat is point1 in form of c(x,y)
#oLonLat is point2 in form of c(x1,y1)
        earthRadius = 6378 #kilometers
        cLon <- cLonLat[1] * (pi/180)
        cLat <- cLonLat[2] * (pi/180)
        oLon <- oLonLat[1] * (pi/180)
        oLat <- oLonLat[2] * (pi/180)
        distance <- sqrt((sin((oLat - cLat)/2))^2 + cos(cLat) *
cos(oLat) * (sin((oLon - cLon)/2))^2)
        distance <- asin(distance) * 6378 * 1000
        distance * 2
}

Dr. Sawada's equal-area projection that allows adjusted origins and
standard latitudes:

albers.mat.f <- function(xy, origin = c(0, 0), stdlns = c(0, 0))
{
#R - Earth Radius
#phi1 - Lat; Stdln 1
#phi2 - Lat; Stdln 2
#phi0 - Lat; Origin
#lambda0 - Lon; Origin
#xy - c(phi,lambda); point to be projected
#C = (cos(phi1))^2+2nsin(phi1)
#n=(sin(phi1)+sin(phi2))/2
#theta = n(lambda-lambda0)
#p = (R(C-2n sin(phi))^.5)/n
#pnot = (R(C-2n sin (phi0))^.5)/n
#x = p sin(theta)
#y = pnot - p cos(theta)
 xy <- as.matrix(xy)
 dimnames(xy) <- list(NULL, c("x", "y")) #dimnames(xy) <- NULL
 xy <- (xy * pi)/180
 R <- 6378 * 1000
 phi1 <- (stdlns[1] * pi)/180
 phi2 <- (stdlns[2] * pi)/180
 phi0 <- (origin[2] * pi)/180
 lambda0 <- (origin[1] * pi)/180 #
############
 n <- (sin(phi1) + sin(phi2))/2
 twon <- 2 * n
 C <- ((cos(phi1))^2 + 2 * n * sin(phi1))
 theta <- n * (xy[, 1] - lambda0)
 xy[, 1] <- ((R * (C - twon * sin(xy[, 2]))^0.5)/n) * sin(theta)
 xy[, 2] <- ((R * (C - twon * sin(phi0))^0.5)/n) - ((R * (C - 2 * n *
sin(xy[, 2]))^0.5)/n) * cos(theta)
 xy
}


--

------------------

Jason Pither
PhD Candidate
c/o Aarssen Lab
Department of Biology
Queen's University
Kingston, ON, CANADA
K7L 3N6
Tel: (h) 613-539-0996
     (w) 613-533-6000 x77327
http://biology.queensu.ca/~pitherj/
pitherj@biology.queensu.ca