Skip to content

LL from bearing and distance

4 messages · Eric Archer, Roger Bivand, Arien Lam

Original
Eric Archer
# 
Is there a package that contains a function that will allow me to 
calculate the ending latitude and longitude given a starting point on 
the earth's surface and bearing and distance?  I have searched the 
archives but have been unable to find one and would like to double check 
before I try to roll my own.  Thanks for any pointers.

Cheers,
eric
Roger Bivand
#
On Thu, 27 Sep 2007, Eric Archer wrote:

            

      
      
      
    

        
Nothing for this case - azimuth for two points is in gzAzimuth() in 
maptools, so that might be a way in. There is C code for the spDistsN1() 
function in the sp package, visible in the CVS repository at the r-spatial 
site at sourceforge, but neither of these meets your need directly. If you 
do "roll your own", please consider contributing to maptools.

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Arien Lam
    

    

      
      #
    

  


  

      
The following function does what Eric describes, if the earth's surface 
is close enough to a sphere.
Roger, I'll try to document it properly for inclusion in Maptools.

Cheers, Arien




# compute the place where you end up, if you travel
# a certain distance along a great circle,
# which is uniquely defined by a point (your starting point)
# and an angle with the meridian at that point (your direction).
# the travelvector is a dataframe with at least the columns
# magnitude and direction.
# n.b. earth radius is the "ellipsoidal quadratic mean radius"
# of the earth, in m.

vectordestination <- function(lonlatpoint, travelvector) {
      Rearth <- 6372795
      Dd <- travelvector$magnitude / Rearth
      Cc <- travelvector$direction

      if (class(lonlatpoint) == "SpatialPoints") {
          lata <- coordinates(lonlatpoint)[1,2] * (pi/180)
          lona <- coordinates(lonlatpoint)[1,1] * (pi/180)
      }
      else {
          lata <- lonlatpoint[2] * (pi/180)
          lona <- lonlatpoint[1] * (pi/180)
      }
      latb <- asin(cos(Cc) * cos(lata) * sin(Dd) + sin(lata) * cos(Dd))
      dlon <- atan2(cos(Dd) - sin(lata) * sin(latb), sin(Cc) * sin(Dd) * 
cos(lata))
      lonb <- lona - dlon + pi/2

      lonb[lonb >  pi] <- lonb[lonb >  pi] - 2 * pi
      lonb[lonb < -pi] <- lonb[lonb < -pi] + 2 * pi

      latb <- latb * (180 / pi)
      lonb <- lonb * (180 / pi)

      cbind(longitude = lonb, latitude = latb)
}

        
Roger Bivand wrote:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Eric Archer
    

    

      
      #
    

  


  

      
All,

I have ported JavaScript code from a couple of websites to R for three 
functions that solve my problem.

'calc.destination.sphere'  assumes a perfect spheroid with a given radius
'calc.destination.ellipsoid' assumes an ellipsoid
'calc.destination.Vincenty' uses the Vincenty model which appears to 
have the most accuracy

I'm unfamiliar with the steps involved in including these functions into 
a package, but will be happy to do it if someone can point me in the 
right direction. 

Cheers,
eric.

---

as.radians <- function(degrees) degrees * pi / 180
as.degrees <- function(radians) radians * 180 / pi
as.bearing <- function(radians) (as.degrees(radians) + 360) %% 360

ellipsoid <- function(model = "WGS84") {
  switch(model,
    WGS84 = c(a = 6378137, b = 6356752.3142, f = 1 / 298.257223563),
    GRS80 = c(a = 6378137, b = 6356752.3141, f = 1 / 298.257222101),
    Airy = c(a = 6377563.396, b = 6356256.909, f = 1 / 299.3249646),
    International = c(a = 6378888, b = 6356911.946, f = 1 / 297),
    Clarke = c(a = 6378249.145, b = 6356514.86955, f = 1 / 293.465),
    GRS67 = c(a = 6378160, b = 6356774.719, f = 1 / 298.25),
    c(a = NA, b = NA, f = NA)
  )
}

calc.destination.sphere <- function(lat, lon, brng, dist, dist.units = 
"km", radius = 6371) {
  # lat, lon : lattitude and longitude in decimal degrees
  # brng : bearing from 0 to 360 degrees
  # dist : distance travelled
  # dist.units : units of distance "km" (kilometers), "nm" (nautical 
miles), "mi" (statute miles)
  # radius : radius of sphere in kilometers
  #
  # Code adapted from JavaScript by Chris Veness 
(scripts at movable-type.co.uk) at 
http://www.movable-type.co.uk/scripts/latlong.html#ellipsoid
  #   originally from Ed Williams' Aviation Formulary, 
http://williams.best.vwh.net/avform.htm#LL
  #
  #

  dist <- switch(dist.units,
    km = dist,
    nm = dist * 1.852,
    mi = dist * 1.609344
  )
  lat <- as.radians(lat)
  lon <- as.radians(lon)
  brng <- as.radians(brng)

  psi <- dist / radius
  lat2 <- asin(sin(lat) * cos(psi) +  cos(lat) * sin(psi) * cos(brng))
  lon2 <- lon + atan2(sin(brng) * sin(psi) * cos(lat), cos(psi) - 
sin(lat) * sin(lat2))
  if (is.nan(lat2) || is.nan(lon2)) return(c(lat = NA, lon = NA))
  c(lat = as.degrees(lat2), lon = as.degrees(lon2))
}

calc.destination.ellipsoid <- function(lat, lon, brng, dist, dist.units 
= "km", model = "WGS84") { 
  # lat, lon : lattitude and longitude in decimal degrees
  # brng : bearing from 0 to 360 degrees
  # dist : distance travelled
  # dist.units : units of distance "km" (kilometers), "nm" (nautical 
miles), "mi" (statute miles)
  # model : choice of ellipsoid model ("WGS84", "GRS80", "Airy", 
"International", "Clarke", "GRS67")
  #
  # Code adapted from JavaScript by Larry Bogan (larry at go.ednet.ns.ca) 
at http://www.go.ednet.ns.ca/~larry/bsc/jslatlng.html
  #
  #

  dist <- switch(dist.units,
    km = dist,
    nm = dist * 1.852,
    mi = dist * 1.609344
  )
  lat <- as.radians(lat)
  lon <- as.radians(lon)
  brng <- as.radians(brng)
  e <- 0.08181922
  ellips <- ellipsoid(model)

  radius <- (ellips["a"] / 1000) * (1 - e ^ 2) / ((1 - e ^ 2 * sin(lat) 
^ 2) ^ 1.5)
  psi <- dist / radius
  phi <- pi / 2 - lat
  arc.cos <- cos(psi) * cos(phi) + sin(psi) * sin(phi) * cos(brng)
  lat2 <- as.degrees((pi / 2) - acos(arc.cos))

  arc.sin <- sin(brng) * sin(psi) / sin(phi)
  lon2 <- as.degrees(lon + asin(arc.sin))
 
  c(lat = lat2, lon = lon2)
}


calc.destination.Vincenty <- function(lat, lon, brng, dist, dist.units = 
"km", model = "WGS84") {
  # lat, lon : lattitude and longitude in decimal degrees
  # brng : bearing from 0 to 360 degrees
  # dist : distance travelled
  # dist.units : units of distance "km" (kilometers), "nm" (nautical 
miles), "mi" (statute miles)
  # model : choice of ellipsoid model ("WGS84", "GRS80", "Airy", 
"International", "Clarke", "GRS67")
  #
  # Code adapted from JavaScript by Chris Veness 
(scripts at movable-type.co.uk) at 
http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html
  # Original reference (http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf):
  #   Vincenty, T. 1975.  Direct and inverse solutions of geodesics on 
the ellipsoid with application of nested equations.
  #      Survey Review 22(176):88-93
  #
  #

  dist <- switch(dist.units,
    km = dist,
    nm = dist * 1.852,
    mi = dist * 1.609344
  )
  dist <- dist * 1000
  lat <- as.radians(lat)
  lon <- as.radians(lon)
  brng <- as.radians(brng)
  ellips <- ellipsoid(model)

  sin.alpha1 <- sin(brng)
  cos.alpha1 <- cos(brng)
  tan.u1 <- (1 - ellips["f"]) * tan(lat)
  cos.u1 <- 1 / sqrt(1 + (tan.u1 ^ 2))
  sin.u1 <- tan.u1 * cos.u1
  sigma1 <- atan2(tan.u1, cos.alpha1)
  sin.alpha <- cos.u1 * sin.alpha1
  cos.sq.alpha <- 1 - (sin.alpha ^ 2)
  u.sq <- cos.sq.alpha * ((ellips["a"] ^ 2) - (ellips["b"] ^ 2)) / 
(ellips["b"] ^ 2)
  cap.A <- 1 + u.sq / 16384 * (4096 + u.sq * (-768 + u.sq * (320 - 175 * 
u.sq)))
  cap.B <- u.sq / 1024 * (256 + u.sq * (-128 + u.sq * (74 - 47 * u.sq)))
 
  sigma <- dist / (ellips["b"] * cap.A)
  sigma.p <- 2 * pi
  cos.2.sigma.m <- cos(2 * sigma1 + sigma)
  while(abs(sigma - sigma.p) > 1e-12) {
    cos.2.sigma.m <- cos(2 * sigma1 + sigma)
    sin.sigma <- sin(sigma)
    cos.sigma <- cos(sigma)
    delta.sigma <- cap.B * sin.sigma * (cos.2.sigma.m + cap.B / 4 * 
(cos.sigma *
       (-1 + 2 * cos.2.sigma.m ^ 2) - cap.B / 6 * cos.2.sigma.m *
       (-3 + 4 * sin.sigma ^ 2) * (-3 + 4 * cos.2.sigma.m ^ 2)))
    sigma.p <- sigma
    sigma <- dist / (ellips["a"] * cap.A) + delta.sigma
  }

  tmp <- sin.u1 * sin.sigma - cos.u1 * cos.sigma * cos.alpha1
  lat2 <- atan2(sin.u1 * cos.sigma + cos.u1 * sin.sigma * cos.alpha1,
    (1 - ellips["f"]) * sqrt(sin.alpha ^ 2 + tmp ^ 2))
  lambda <- atan2(sin.sigma * sin.alpha1, cos.u1 * cos.sigma - sin.u1 * 
sin.sigma * cos.alpha1)
  cap.C <- ellips["f"] / 16 * cos.sq.alpha * (4 + ellips["f"] * 
(ellips["f"] - 3 * cos.sq.alpha))
  cap.L <- lambda - (1 - cap.C) * ellips["f"] * sin.alpha *
    (sigma + cap.C * sin.sigma * (cos.2.sigma.m + cap.C * cos.sigma * 
(-1 + 2 * cos.2.sigma.m ^ 2)))
 
  c(lat = as.degrees(lat2), lon = as.degrees(lon + cap.L))
}