function to generate new location given a lon/lat bearing and distance

get_point_at_distance.R

get_point_at_distance <- function(lon, lat, d, bearing, R = 6378137) {
  # lat: initial latitude, in degrees
  # lon: initial longitude, in degrees
  # d: target distance from initial point (in m)
  # bearing: (true) heading in degrees
  # R: mean radius of earth (in m)
  # Returns new lat/lon coordinate {d} m from initial, in degrees
  stopifnot("lon value not in range {-180, 180}" = !any(lon < -180 | lon > 180))
  stopifnot("lat value not in range {-90, 90}" = !any(lat < -90 | lat > 90))
  stopifnot("bearing not in range {0, 360}" = !any(bearing < 0 | bearing > 360))
  warning("I recommend doing this in projected coordinates at small scales!",
          immediate. = TRUE, call. = FALSE)
  ## convert to radians
  lat1 <- lat * (pi/180)
  lon1 <- lon * (pi/180)
  a <- bearing * (pi/180)
  ## new position
  lat2 <- asin(sin(lat1) * cos(d/R) + cos(lat1) * sin(d/R) * cos(a))
  lon2 <- lon1 + atan2(
    sin(a) * sin(d/R) * cos(lat1),
    cos(d/R) - sin(lat1) * sin(lat2)
  )
  ## convert back to degrees
  lat2 <- lat2 * (180/pi)
  lon2 <- lon2 * (180/pi)
  ## return
  return(cbind(lon2, lat2))
}

TEST <- FALSE
if(TEST) {
  lat = 52.20472
  lon = 0.14056
  distance = 15000
  bearing = 90
  get_point_at_distance(lon = lon, lat = lat, d = distance, bearing = bearing)
  # [1]  0.3604322 52.2045157
}


get_point_at_distance_projected <- function(X, Y, d, bearing) {
  # X: projected X coordinate
  # Y: projected Y coordinate
  # d: target distance from initial point (in m)
  # bearing: (true) heading in degrees
  # Returns new X/Y coordinate {d} m from initial, in degrees
  stopifnot("bearing not in range {0, 360}" = !any(bearing < 0 | bearing > 360))
  bearing_rad <- bearing * (pi / 180)
  delta_x <- d * sin(bearing_rad)
  delta_y <- d * cos(bearing_rad)
  new_x <- X + delta_x
  new_y <- Y + delta_y
  return(cbind(new_x, new_y))
}

if(TEST) {
  X = 331137
  Y = 6158080
  d = 745
  bearing = 90
  get_point_at_distance_projected(X, Y, d, bearing)
}