Jason:
Here is my great circle distance function in S+
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
}
I also use Albers Equal Area Conic projection often (though not
equidistant it is allright for my purposes) and here is the function
which I wrote to overcome the fact that the map projections built into
Splus don't allow chaning origins or standard lines:
albers.mat.f<-function(xy, origin = c(0, 0), stdlns = c(0, 0)) {
#inputs############# #xy - a two column dataframe or matrix with
longitudes in the first column (x) and latitudes in the second column
#(y) #origin = c(0,0): is the Longitude and Latitude for the center of
the output coordinate system e.g., c(-100,50)
# for -100 West and 50 N.
#stdlns = c(0,0) are the standard lines (parallels) where true scale is
found e.g., for North America perhaps
#c(35,80) for 35 N and 80 N since this distributes distortion throughout
the map. #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
}
You can send your data to a function like this before sending to the
neighbor function of splus spatial stasts. Other equations are easily
computed from standard projection equations.
regards,
Dr. M.Sawada
Assistant Professor, GIS
University of Ottawa
Dept. Geography
P.O. Box 450 Stn. A
Ottawa, ON K1N 6N9
Tel: 613-562-5800 x1040
Fax: 613-562-5145
Eml: msawada@uottawa.ca
-----Original Message-----
From: s-news-owner@lists.biostat.wustl.edu
[mailto:s-news-owner@lists.biostat.wustl.edu] On Behalf Of Jason Pither
Sent: Tuesday, December 11, 2001 1:52 PM
To: s-news@lists.biostat.wustl.edu
Subject: [S] great circle distance and find.neighbor
Hello:
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.
Any suggestions would be much appreciated.
Jason Pither.
--
------------------
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
---------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu. To
unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
the BODY of the message: unsubscribe s-news
|