[R-sig-Geo] function for bearing
Roger Bivand
Roger.Bivand at nhh.no
Thu Mar 9 18:06:48 CET 2006
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: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list