Plotting GPS tracks with R

Introduction

This is a notebook based on a blog written by Sascha Wolfer.

Many GPS devices and apps have the capability to track your current position via GPS. If you go walking, running, cycling, flying or driving, you can take a look at your exact route and your average speed. Some of these devices or apps also allow you to export your routes in various formats, e.g., the popular XML-based GPX format.

The Sascha’s notebook shows how to: – read in a GPX file using R and its XML package, – calculate distances and speeds between points, – plot elevation and speed, – plot a track, – plot a track on a map.

Here are some of Sacha’s results:

• The altitude of the track in time

• The speed during the run, with smoother

• The raw track without a map in the background

• The track with an OpenStreetMap of Stuttgart in the background

• The track with a Bing map of Stuttgart in the background

We will replicate some of the code (and will add new code if issues arise).

The code

First, we load several libraries:

library(XML)
library(OpenStreetMap)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:base':
## 
##     date

library(ggmap)

## Loading required package: ggplot2

library(ggplot2)
library(raster)

## Loading required package: sp

library(sp)

Now, we define a function that shifts vectors conveniently:

shift.vec <- function (vec, shift) {
  if(length(vec) <= abs(shift)) {
    rep(NA ,length(vec))
  }else{
    if (shift >= 0) {
      c(rep(NA, shift), vec[1:(length(vec)-shift)]) }
    else {
      c(vec[(abs(shift)+1):length(vec)], rep(NA, abs(shift))) } } }

Before going, let’s check what this function is for:

col1 <- seq(0,100,5)

col2 <- seq(200, 100, -5)

my_df <- data.frame(c1= col1, c2= col2)

my_df

my_df$nc1 <- shift.vec(my_df$c1, -1)
my_df$nc2 <- shift.vec(my_df$c2, -1)
my_df

Now, we’re reading in the GPX file that Sacha shared here. After downloading the GPX file, we have to upload it to our RStudio Cloud project. Then, we can proceed:

options(digits=10)
# Parse the GPX file
pfile <- htmlTreeParse(file = "./datos/run2.gpx", error = function(...) {
}, useInternalNodes = T)

# Get all elevations, times and coordinates via the respective xpath
elevations <- as.numeric(xpathSApply(pfile, path = "//trkpt/ele", xmlValue))
times <- xpathSApply(pfile, path = "//trkpt/time", xmlValue)
coords <- xpathSApply(pfile, path = "//trkpt", xmlAttrs)

str(coords)

##  chr [1:2, 1:931] "48.78456666666666" "9.2173" "48.78456666666666" ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:2] "lat" "lon"
##   ..$ : NULL

# Extract latitude and longitude from the coordinates
lats <- as.numeric(coords["lat",])
lons <- as.numeric(coords["lon",])

# Put everything in a dataframe and get rid of old variables
geodf <- data.frame(lat = lats, lon = lons, ele = elevations, time = times)
rm(list=c("elevations", "lats", "lons", "pfile", "times", "coords"))
head(geodf)

We already have nice dataframe now with all the information available in the GPX file for each position. Each position is defined by the latitude and longitude and we also have the elevation (altitude) available. Note, that the altitude ist quite noisy with GPS, we will see this in a minute.

Now, let’s calculate the distances between successive positions and the respective speed in this segment.

# Shift vectors for lat and lon so that each row also contains the next position.
geodf$lat.p1 <- shift.vec(geodf$lat, -1)
geodf$lon.p1 <- shift.vec(geodf$lon, -1)
head(geodf)

# Calculate distances (in metres) using the function pointDistance from the ‘raster’ package.
# Parameter ‘lonlat’ has to be TRUE!

geodf$dist.to.prev <- apply(geodf, 1, FUN = function (row) {
  pointDistance(c(as.numeric(row["lat.p1"]),
  as.numeric(row["lon.p1"])),
                c(as.numeric(row["lat"]), as.numeric(row["lon"])),
                lonlat = T)
})

head(geodf$dist.to.prev)

## [1]  0.00000000 10.27107832 10.33552017 10.05242493 16.75869953 14.36471901

td <- sum(geodf$dist.to.prev, na.rm=TRUE)
print(paste("The distance run was ", td, " meters"))

## [1] "The distance run was  11097.0963172583  meters"

# Transform the column ‘time’ so that R knows how to interpret it.
geodf$time <- strptime(geodf$time, format = "%Y-%m-%dT%H:%M:%OS")
# Shift the time vector, too.
geodf$time.p1 <- shift.vec(geodf$time, -1)
# Calculate the number of seconds between two positions.
geodf$time.diff.to.prev <- as.numeric(difftime(geodf$time.p1, geodf$time))

head(geodf$time.diff.to.prev, n=15) 

##  [1] 0.06999993324 3.00999999046 2.97000002861 6.00999999046 2.98000001907
##  [6] 2.09999990463 1.93000006676 1.96000003815 2.05999994278 3.01999998093
## [11] 2.94000005722 3.07999992371 2.94000005722 3.87000012398 2.19999980927

# Calculate metres per seconds, kilometres per hour and two LOWESS smoothers to get rid of some noise.
geodf$speed.m.per.sec <- geodf$dist.to.prev / geodf$time.diff.to.prev
geodf$speed.km.per.h <- geodf$speed.m.per.sec * 3.6
geodf$speed.km.per.h <- ifelse(is.na(geodf$speed.km.per.h), 0, geodf$speed.km.per.h)
geodf$lowess.speed <- lowess(geodf$speed.km.per.h, f = 0.2)$y
geodf$lowess.ele <- lowess(geodf$ele, f = 0.2)$y

Now, let’s plot all the stuff!

# Plot elevations and smoother
plot(geodf$ele, type = "l", bty = "n", xaxt = "n", ylab = "Elevation", xlab = "", col = "grey40")
lines(geodf$lowess.ele, col = "red", lwd = 3)
legend(x="bottomright", legend = c("GPS elevation", "LOWESS elevation"),
       col = c("grey40", "red"), lwd = c(1,3), bty = "n")

# Plot speeds and smoother
plot(geodf$speed.km.per.h, type = "l", bty = "n", xaxt = "n", ylab = "Speed (km/h)", xlab = "",
     col = "grey40")
lines(geodf$lowess.speed, col = "blue", lwd = 3)
legend(x="bottom", legend = c("GPS speed", "LOWESS speed"),
       col = c("grey40", "blue"), lwd = c(1,3), bty = "n")
abline(h = mean(geodf$speed.km.per.h), lty = 2, col = "blue")

# Plot the track without any map, the shape of the track is already visible.
plot(rev(geodf$lon), rev(geodf$lat), type = "l", col = "red", lwd = 3, bty = "n", ylab = "Latitude", xlab = "Longitude")

We will now use the OpenStreetMap package to get some maps from the internet and use them as a background for the track we just converted. There are several blocks for each map type. Check the comments in the first block what the different calls do.

First, we need to get the specific map. The type argument controls which type of map we get.

Now, let´s plot the track using the ggmap package:

library(ggmap)
lat <- c(min(geodf$lat), max(geodf$lat))
lat

## [1] 48.78401500 48.80601333

lon <- c(min(geodf$lon), max(geodf$lon))
lon

## [1] 9.191535 9.221930

bbox <- make_bbox(lon,lat)

b1<- get_map(bbox,maptype="watercolor", source="stamen")

## Map from URL : http://tile.stamen.com/watercolor/15/17220/11279.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17221/11279.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17222/11279.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17223/11279.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17220/11280.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17221/11280.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17222/11280.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17223/11280.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17220/11281.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17221/11281.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17222/11281.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17223/11281.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17220/11282.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17221/11282.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17222/11282.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17223/11282.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17220/11283.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17221/11283.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17222/11283.jpg

## Map from URL : http://tile.stamen.com/watercolor/15/17223/11283.jpg

ggmap(b1) + geom_point(data = geodf, 
           aes(lon,lat,col = ele), size=1, alpha=0.7) +
           labs(x = "Longitude", y = "Latitude",
           title="Track of one Sascha run through Stuttgart")

Another map, now using the mapview package. It was created to help researchers during their spatial data analysis workflow. It provides functions to very quickly and conveniently create interactive visualisations of spatial data.

For such a task, we need to convert the data frame into a spatial data frame. The functions coordinates and projection from sp library converts the dataset into spatial objects that mapview supports.

In particular, coordinates specifies the latitude and longitude of the data, and proj4string creates the projection layer, i.e. the coordinate system. As you know, the EPSG: 4326 means that the coordinates are latitude-longitude pairs on a reference ellipsoid given by WGS84, projected with Mercator Projection.

library(mapview)

class(geodf)

## [1] "data.frame"

spdf_geo <- geodf

coordinates(spdf_geo) <- ~ lon + lat
proj4string(spdf_geo) <- "+init=epsg:4326"

class(spdf_geo)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

mapview(spdf_geo)

Finally, we will use leaflet for mapping the track:

library(leaflet)

leaflet() %>% 
  addTiles() %>% 
  addFeatures(spdf_geo, weight = 1, fillColor = "grey", color = "black",
              opacity = 1, fillOpacity = 0.6)

That’s all for now.