On Thu, 9 Mar 2006, Sebastian Luque wrote: > Hello, > > I have searched the archives and other sources, but cannot find an R > function that would calculate the bearing of a pair of points (lat, lon). > It wouldn't be difficult to write such a function, but has this already > been done in a package?
I've an implementation that I will add to sp or maptools. It calculates the azimuth using Snyder's spherical equations (p. 30), or an equivalent solution by Abdali: gzAzimuth <- function(from, to, type="snyder_sphere") { deg2rad <- function(x) x*pi/180 rad2deg <- function(x) x*180/pi # note negative longitudes if (is.matrix(from)) { lon <- -deg2rad(from[,1]) lat <- deg2rad(from[,2]) } else { lon <- -deg2rad(from[1]) lat <- deg2rad(from[2]) } if (is.matrix(to)) { if (nrow(to) > 1) stop("to: single coordinate") to <- c(to) } else { lon0 <- -deg2rad(to[1]) lat0 <- deg2rad(to[2]) } dflon = lon-lon0 # results in degrees from N, negative west if (type == "abdali") res <- atan2(sin(dflon), ((cos(lat)*tan(lat0)) - (sin(lat)*cos(dflon)))) else if (type == "snyder_sphere") res <- atan2((cos(lat0)*sin(dflon)), (cos(lat)*sin(lat0)) - (sin(lat)*cos(lat0)*cos(dflon))) else stop("type unkown") is.na(res) <- lon == lon0 & lat == lat0 rad2deg(res) } There is a very interesting discussion of the centrality of azimuth-finding in the development of mathematics and mathematical geography in: http://patriot.net/users/abdali/ftp/qibla.pdf. Among others, al-Khwarizmi was an important contributor. It seems to work like Tim Keitt's suggestion: > gzAzimuth(c(10,50), c(70,20)) [1] 99.76647 > 99 + (45 + (59.287/60))/60 [1] 99.76647 $ geod +ellps=sphere << EOF -I > 50N 10E 20N 70E > EOF 99d45'59.287" -42d23'10.939" 6189790.454 The first argument to the function can be a matrix of coordinates, the second is a single coordinate. Roger > > > Thanks in advance, > > -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [EMAIL PROTECTED] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo