1# 2007 Eric Archer and Roger Bivand
2#
3
4gcDestination <- function(lon, lat, bearing, dist, dist.units = "km",
5    model=NULL, Vincenty=FALSE) {
6 # lat, lon : lattitude and longitude in decimal degrees
7 # bearing : bearing from 0 to 360 degrees
8 # dist : distance travelled
9 # dist.units : units of distance "km" (kilometers), "nm" (nautical
10 # miles), "mi" (statute miles)
11 # model : choice of ellipsoid model ("WGS84", "GRS80", "Airy",
12 # "International", "Clarke", "GRS67")
13
14    if (!is.numeric(lon)) stop("lon not numeric")
15    if (!is.numeric(lat)) stop("lat not numeric")
16    if (!is.numeric(bearing)) stop("bearing not numeric")
17    if (!is.numeric(dist)) stop("dist not numeric")
18
19    if (length(lon) != length(lat)) stop("lon and lat differ in length")
20    if (length(bearing) > 1L && length(lon) > 1L) stop("length mismatch")
21    if (length(bearing) > 1L && length(dist) > 1L) stop("length mismatch")
22
23    as.radians <- function(degrees) degrees * pi / 180
24    as.degrees <- function(radians) radians * 180 / pi
25    as.bearing <- function(radians) (as.degrees(radians) + 360) %% 360
26
27    ellipsoid <- function(model = "WGS84") {
28        switch(model,
29        WGS84 = c(a = 6378137, b = 6356752.3142, f = 1 / 298.257223563),
30        GRS80 = c(a = 6378137, b = 6356752.3141, f = 1 / 298.257222101),
31        Airy = c(a = 6377563.396, b = 6356256.909, f = 1 / 299.3249646),
32        International = c(a = 6378888, b = 6356911.946, f = 1 / 297),
33        Clarke = c(a = 6378249.145, b = 6356514.86955, f = 1 / 293.465),
34        GRS67 = c(a = 6378160, b = 6356774.719, f = 1 / 298.25),
35        c(a = NA, b = NA, f = NA)
36    )}
37
38    dist <- switch(dist.units,
39        km = dist,
40        nm = dist * 1.852,
41        mi = dist * 1.609344
42    )
43    lat <- as.radians(lat)
44    lon <- as.radians(lon)
45    bearing <- as.radians(bearing)
46
47    if (is.null(model)) {
48 # Code adapted from JavaScript by Chris Veness
49 # (scripts@movable-type.co.uk) at
50 # http://www.movable-type.co.uk/scripts/latlong.html#ellipsoid
51 #   originally from Ed Williams' Aviation Formulary,
52 # http://williams.best.vwh.net/avform.htm
53        radius <- 6371
54        psi <- dist / radius
55        lat2 <- asin(sin(lat) * cos(psi) +  cos(lat) * sin(psi) * cos(bearing))
56        lon2 <- lon + atan2(sin(bearing) * sin(psi) * cos(lat), cos(psi) -
57            sin(lat) * sin(lat2))
58        if (any(is.nan(lat2)) || any(is.nan(lon2))) warning("Out of range values")
59        return(cbind(long=as.degrees(lon2), lat=as.degrees(lat2)))
60    }
61
62    ellips <- ellipsoid(model)
63    if (is.na(ellips["a"])) stop("no such ellipsoid model")
64    if (Vincenty) {
65 # Code adapted from JavaScript by Chris Veness
66 # (scripts@movable-type.co.uk) at
67 # http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html
68 # Original reference (http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf):
69 #   Vincenty, T. 1975.  Direct and inverse solutions of geodesics on
70 # the ellipsoid with application of nested equations.
71 #      Survey Review 22(176):88-93
72        dist <- dist * 1000
73        sin.alpha1 <- sin(bearing)
74        cos.alpha1 <- cos(bearing)
75        tan.u1 <- (1 - ellips["f"]) * tan(lat)
76        cos.u1 <- 1 / sqrt(1 + (tan.u1 ^ 2))
77        sin.u1 <- tan.u1 * cos.u1
78        sigma1 <- atan2(tan.u1, cos.alpha1)
79        sin.alpha <- cos.u1 * sin.alpha1
80        cos.sq.alpha <- 1 - (sin.alpha ^ 2)
81        u.sq <- cos.sq.alpha * ((ellips["a"] ^ 2) - (ellips["b"] ^ 2)) /
82            (ellips["b"] ^ 2)
83        cap.A <- 1 + u.sq / 16384 * (4096 + u.sq * (-768 + u.sq * (320 -
84            175 * u.sq)))
85        cap.B <- u.sq / 1024 * (256 + u.sq * (-128 + u.sq * (74 - 47 * u.sq)))
86
87        sigma <- dist / (ellips["b"] * cap.A)
88        sigma.p <- 2 * pi
89        cos.2.sigma.m <- cos(2 * sigma1 + sigma)
90        while(any(abs(sigma - sigma.p) > 1e-12)) {
91            cos.2.sigma.m <- cos(2 * sigma1 + sigma)
92            sin.sigma <- sin(sigma)
93            cos.sigma <- cos(sigma)
94            delta.sigma <- cap.B * sin.sigma * (cos.2.sigma.m + cap.B / 4 *
95                (cos.sigma *
96                (-1 + 2 * cos.2.sigma.m ^ 2) - cap.B / 6 * cos.2.sigma.m *
97                (-3 + 4 * sin.sigma ^ 2) * (-3 + 4 * cos.2.sigma.m ^ 2)))
98            sigma.p <- sigma
99            sigma <- dist / (ellips["a"] * cap.A) + delta.sigma
100        }
101        tmp <- sin.u1 * sin.sigma - cos.u1 * cos.sigma * cos.alpha1
102        lat2 <- atan2(sin.u1 * cos.sigma + cos.u1 * sin.sigma * cos.alpha1,
103            (1 - ellips["f"]) * sqrt(sin.alpha ^ 2 + tmp ^ 2))
104        lambda <- atan2(sin.sigma * sin.alpha1, cos.u1 * cos.sigma - sin.u1 *
105            sin.sigma * cos.alpha1)
106        cap.C <- ellips["f"] / 16 * cos.sq.alpha * (4 + ellips["f"] *
107            (ellips["f"] - 3 * cos.sq.alpha))
108        cap.L <- lambda - (1 - cap.C) * ellips["f"] * sin.alpha *
109            (sigma + cap.C * sin.sigma * (cos.2.sigma.m + cap.C * cos.sigma *
110            (-1 + 2 * cos.2.sigma.m ^ 2)))
111        lat2 <- as.degrees(lat2)
112        lon2 <- as.degrees(lon + cap.L)
113    } else {
114 # Code adapted from JavaScript by Larry Bogan (larry@go.ednet.ns.ca)
115 # at http://www.go.ednet.ns.ca/~larry/bsc/jslatlng.html
116        e <- 0.08181922
117        radius <- (ellips["a"] / 1000) * (1 - e^2) / ((1 - e^2 *
118            sin(lat)^2)^1.5)
119        psi <- dist / radius
120        phi <- pi / 2 - lat
121        arc.cos <- cos(psi) * cos(phi) + sin(psi) * sin(phi) * cos(bearing)
122        lat2 <- as.degrees((pi / 2) - acos(arc.cos))
123        arc.sin <- sin(bearing) * sin(psi) / sin(phi)
124        lon2 <- as.degrees(lon + asin(arc.sin))
125    }
126    return(cbind(long=lon2, lat=lat2))
127}
128