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