[R-sig-Geo] LL from bearing and distance

Arien Lam a.lam at geo.uu.nl
Thu Sep 27 23:10:21 CEST 2007


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:
> On Thu, 27 Sep 2007, Eric Archer wrote:
> 
>> 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.
> 
> 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
> 
>> Cheers,
>> eric
>>
>>
>




More information about the R-sig-Geo mailing list