I want to get the coordinates of new points, preferably using the sf package, when the inital positions and the distance and heading traveled is known. Consider this; we have three points (pts), with a heading and a distance in km attached. How to find the coordinates for the new positions?
library(data.table, dplyr, sf)
dat <- data.table(lon = c(10,10.1,10.4), lat = c(58.4,57.4,57.8),
heading = c(45,10,235), distance_km = c(1,5.3,3))
pts <- dat %>%
sf::st_as_sf(coords = c("lon","lat")) %>%
sf::st_set_crs(4326)
Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 10 ymin: 57.4 xmax: 10.4 ymax: 58.4
Geodetic CRS: WGS 84
heading distance_km geometry
1 45 1.0 POINT (10 58.4)
2 10 5.3 POINT (10.1 57.4)
3 235 3.0 POINT (10.4 57.8)
Was considering making circles around the points, but dont know how to connect the point to the circle with the correct heading.
buf <- st_buffer(pts, dist = pts$distance_km*1000)
circ <- st_cast(buf, "LINESTRING")
Found the answer here: Calculate coordinate from starting point, having distance and an angle for all quadrants and here: Convert radians to degree / degree to radians
Will post my R solution for completeness. If anyone has a better or more smooth solution, feel free to post it.
library(data.table, sf, mapview)
dat <- data.table(lon = c(10,10.1,10.4), lat = c(58.4,57.4,57.8),
heading = c(45,10,235), distance_km = c(1,5.3,3))
pts <- dat %>%
sf::st_as_sf(coords = c("lon","lat")) %>%
sf::st_set_crs(4326)
pts <- st_transform(pts, 32632)
pts$utm_n <- st_coordinates(pts)[,1]
pts$utm_e <- st_coordinates(pts)[,2]
buf <- st_buffer(pts, dist = pts$distance_km*1000)
circ <- st_cast(buf, "LINESTRING")
rad2deg <- function(rad) {(rad * 180) / (pi)}
deg2rad <- function(deg) {(deg * pi) / (180)}
pts$newp_e <- pts$utm_e + (pts$distance_km*1000* cos(deg2rad(pts$heading)))
pts$newp_n <- pts$utm_n + (pts$distance_km*1000* sin(deg2rad(pts$heading)))
dt <- data.table(pts)
pts2 <- dt %>%
sf::st_as_sf(coords = c("newp_n", "newp_e")) %>%
sf::st_set_crs(32632)
mapview(pts2) + mapview(pts, zcol = "heading") + mapview(circ)