0

Using tmaptools package in R - How can I extract the 'Bearing' information from a .GPX track file. This appears in Garmin Basecamp but does not appear using tmaptools::read_GPX. Currently I use the below code. But surely there is a simpler way? Link to GPS Track: https://www.dropbox.com/s/02p3yyjkv9fmrni/Barron_Thomatis_2019_EOD.gpx?dl=0

library(tmaptools)
library(tmap)
library(sf)
library(tidyverse)
library(geosphere)


GPSTrack <- read_GPX("Barron_Thomatis_2019_EOD.gpx", layers = "track_points", as.sf = TRUE)

#
#Adjust GPS Track Data
#

#Extract Lat & Lon from Track geometery (c(lat, Lon))
GPSTrack_Pts <- st_coordinates(GPSTrack)

#Add X, Y Columns to Track
GPSTrack2 <- cbind(GPSTrack, GPSTrack_Pts)

#Create a coordinate vector by combining X & Y
coords <- cbind(GPSTrack2$X,GPSTrack2$Y)

#Convert GPS Track into SpatialPoints format for calculating Bearing
GPSTrack_SpPts <- SpatialPoints(coords)

#Create GPS Point Bearing, GPP point distance & GPS Time interval columns
empty <- st_as_sfc("POINT(EMPTY)")

GPSTrack2 <- GPSTrack2 %>%
  st_set_crs(4326) %>%  # will use great circle distance
  mutate(
    Bearing = bearing(coords))

#Convert Bearing to Course and Add as column
GPSTrack2 <- GPSTrack2 %>% 
  mutate(course = (Bearing + 360) %% 360) # add full circle, i.e. +360, and determine modulo for 360
Jock
  • 75
  • 1
  • 9

1 Answers1

2

I suggest you use lwgeom::st_geod_azimuth() for this task - it makes for somewhat more concise code.

Note that there is a challenge when adding the vector of bearings back to the spatial dataframe of points; it has by definition one element less than is the number of rows (you need two points to define a bearing).

One possibility of achieving that - if required - is by concatenating the vector with a single NA value representing the bearing of the very last point. By definition it has no azimuth, as there is no following point.

The azimuth values are objects of class units, originally in radians. Should the class create a problem (as it does with concatenating with the NA) you can easily convert it to a plain number via units::drop_units().

library(sf)
library(dplyr)
library(lwgeom)


points <- st_read("Barron_Thomatis_2019_EOD.gpx",
                  layer = "track_points",
                  quiet = T,
                  stringsAsFactors = F)

points <- points %>% 
  mutate(bearing = c(lwgeom::st_geod_azimuth(.) %>% units::drop_units(), NA))
Jindra Lacko
  • 7,814
  • 3
  • 22
  • 44
  • Thanks @Jindra Lacko. As I do wish to add this information back to the points spatial dataframe. How can I add the bearing information Point 1 to Point 2 etc with the last point being blank having no bearing? – Jock Mar 06 '20 at 00:09
  • 1
    No problem @Jock! I have edited the code slightly. It is no longer as concise, but the bearing is now a column of the original data frame. With azimuth of the last point set to NA, as there is no next point (the last point being last by definition :) – Jindra Lacko Mar 06 '20 at 07:47