Last active
June 28, 2023 23:58
-
-
Save scbrown86/88a8bfcfcd75e542d3e771b65192aa0f to your computer and use it in GitHub Desktop.
function to generate new location given a lon/lat bearing and distance
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment