Скрещиваем R и велоспорт (анализ спортивных треков)

Я люблю бегать и кататься на велосипеде, и часто возникает желание подумать как проходила та или иная тренировка. Обычно я записываю GPS-трек и просматриваю его потом в сервисах Endomondo, Strava, в которых имеются довольно развитые средства анализа, но раз уж я владею R, то интересно и исследовать треки в этой системе. Ниже я опишу по шагам как это делается.

.1. Библиотеки. Помимо стандартных библиотек для работы с данными и графикой, потребуются библиотеки для загрузки данных в формате XML, в котором сохраняют треки все спортивные устройства и программы вроде Endomondo и Strava, пакеты для расчета расстояний на базе геоданных (geosphere) и рисования при помощи карт Google Map.

###
# install.packages("data.table")
# install.packages("lubridate")
# install.packages("XML")
# install.packages("ggmap")
# install.packages("Imap")
# install.packages("raster")
# install.packages("geosphere")

library(data.table)
library(lubridate)
library(ggplot2)
library(xml2)

library(geosphere)
library(raster)
library(Imap)
library(ggmap)

.2. Для загрузки данных я создал две функции. Первая читает данные в формате TCX (Training Center XML, используется, например, в программе Endomondo), вторая — в формате GPX (Garmin, Strava).

Обе функции интенсивно используют пакет XML2 для разборки XML-файла. Этот пакет использует синтаксис языка XPath для нахождения нужных узлов графа, и извлечения из них информации. Если строка XPath начинается с '\', то эти функции производят глобальный поиск во всем документе, если же нет, то только поиск в пределах уже выбранного узла. Функция xml_find_all выдает список всех узлов, отвечающих критериям, а xml_find_first находит только первый узел. Функция xml_children получает все дочерние узлы выбранного узла, xml_text извлекает из узла его текст, xml_attr — значения атрибутов узла.

Извлеченные данные складываются в список элементов типа data.frame, которые потом быстро объединяются в один data.frame при помощи функции rbindlist.

##Load libraries

load_tcx <- function(fname, gethr = FALSE ) {
  dat1 <- read_xml(fname, options=NULL)

  xml_ns_strip( dat1 )

  #Find the first track subtree
  track1 <- xml_find_first(dat1, "//Track")

#  lapply( xml_find_all(track1, "//Trackpoint"), print )

  tpt <- xml_find_all(track1, "Trackpoint")
  df1 <- rbindlist ( lapply( tpt, function(pt) 
    { 
      return( 
      data.frame( lat = as.numeric( xml_text(xml_find_first(pt, "Position/LatitudeDegrees"))),
        lon = as.numeric( xml_text( xml_find_first(pt, "Position/LongitudeDegrees"))),
        alt = as.numeric( xml_text( xml_find_first(pt, "AltitudeMeters"))),
        t = strptime(xml_text( xml_find_first(pt, "Time")), format = "%Y-%m-%dT%H:%M:%S"),
        hr = as.numeric( xml_text( xml_find_first(pt,"HeartRateBpm/Value"))),
        ds0_m = as.numeric( xml_text( xml_find_first(pt, "DistanceMeters"))),  
        stringsAsFactors = FALSE))})  )   

  return(df1)  
}


load_gpx <- function(fname) {
  dat1 <- read_xml(fname, options=NULL)

  xml_ns_strip( dat1 )

  track1 <- xml_find_first(dat1, "//trkseg")

  df1 <- rbindlist ( lapply(xml_children(track1), function(pt) 
  { return( 
    data.frame( lat = as.numeric( xml_attr( pt, "lat") ),
                lon = as.numeric( xml_attr( pt, "lon") ),
                alt = as.numeric( xml_text( xml_find_first(pt, "ele"))),
                t = strptime(xml_text( xml_find_first(pt, "time")), format = "%Y-%m-%dT%H:%M:%S"),
                hr = as.numeric( xml_text( xml_find_first(pt,"extensions/gpxtpx:TrackPointExtension/gpxtpx:hr"))),
                stringsAsFactors = FALSE))})  )    

  return(df1)  
}

.3. После того как вся информация считана в одну большую таблицу, ее можно расширить вычисляемыми полями. Вычисление скоростей по геокоординатам производится при помощи функции пакета geosphere. Этот пакет предлагает несколько разных методов для подобного расчета, я использовал “Эллипсоид Винсенти” (distVincentyEllipsoid).

Так как данные получаются зашумленными, то для визуализации их стоит сглаживать. Для сглаживания я вначале использовал метод lowess, но эксперименты показали, что в лесу навигаторы иногда измеряют координаты с таким уровнем шума, что lowess дает бессмысленные результаты. Потому я в итоге просто остановился на вычислении средней скорости исходя из дистанции, соответствующей какому-то относительно большому числу точек (параметр nsmooth). Относительно хорошие результаты получаются если это число соответствует промежутку времени ~90 сек.

#Вспомогательные функции
fwd_vec <- function( v ) {
  return( c( v[ -1 ], v[ length(v)] ))
}

lag_vec <- function( v ) {
  return( c( v[1], v[ 1:(length(v)-1)]))
}

#Вычисление средней скорости на отрезке
av_speed <- function( cumdist, cumtime, nsmooth ) {
  v <- diff(cumdist, lag=nsmooth)/diff(cumtime, lag=nsmooth)
  return( c( rep( v[1], floor(nsmooth/2)), v, rep(v[length(v)],nsmooth - floor(nsmooth/2))))
}

#Расчет оптимального числа точек для сглаживания
opt_nsmooth <- function(df){
  #Optimal number of point to smooth - 90 sec 
  return (as.integer(90/as.numeric(mean( diff( df$t)))))

}

#Вычисляемые поля 
calc_spd <- function( df1, nsmooth =50 ) {

  df2 <- df1

  # Time step
  df2$dt <- as.numeric(difftime(df2$t, lag_vec( df2$t)) )

  #Previous points
  df2$lat1 <- lag_vec(df2$lat)
  df2$lon1 <- lag_vec(df2$lon)
  df2$alt1 <- lag_vec(df2$alt)

  # Calculate distances (in metres) using function from the raster package.

  #some woodoo with as.numeric is required  
  df2$ds <- apply(df2, 1, FUN = function (row) {
    distVincentyEllipsoid(c(as.numeric(row[["lon1"]]), as.numeric(row[["lat1"]])),
                          c(as.numeric(row[["lon"]]), as.numeric(row[["lat"]])))
  }) 

# Аналог
  # df2$ds_geo <- apply(df2, 1, FUN = function (row) {
  #   distGeo(c(as.numeric(row[["lon1"]]), as.numeric(row[["lat1"]])),
  #                         c(as.numeric(row[["lon"]]), as.numeric(row[["lat"]])))
  # }) 

#correct distances on elevations

  df2$ds_alt <- sqrt( df2$ds*df2$ds + (df2$alt-df2$alt1)*(df2$alt-df2$alt1))

  df2$ds  <- ifelse(is.na(df2$ds ), 0, df2$ds)
  df2$ds_alt  <- ifelse(is.na(df2$ds_alt ), 0, df2$ds_alt)

  # Calculate speed kilometres per hour and two LOWESS smoothers to get rid of some noise.
  df2$speed_kmh  <- df2$ds / df2$dt  * 3.6
  df2$speed_kmh  <- ifelse(is.na(df2$speed_kmh ), 0, df2$speed_kmh)

  df2$speed_kmh_alt  <- df2$ds_alt / df2$dt  * 3.6
  df2$speed_kmh_alt  <- ifelse(is.na(df2$speed_kmh_alt ), 0, df2$speed_kmh_alt)

  #Smoothed speed
  #Unfortunately lowees sometimes go crazy with GPS data
  #df2$speed_lowess <- lowess(df2$speed_kmh, f = 100/nsmooth)$y
  #df2$speed_lowess_alt <- lowess(df2$speed_kmh_alt, f = 100/nsmooth)$y

  #Distance along the track
  df2$distance <- cumsum( df2$ds)/1000
  df2$distance_alt <- cumsum( df2$ds_alt)/1000

  df2$cumtime <- cumsum( df2$dt )/3600 #hours

  df2$speed_lowess <- av_speed(df2$distance, df2$cumtime, nsmooth)
  df2$speed_lowess_alt <- av_speed(df2$distance_alt, df2$cumtime, nsmooth)

  df2$alt_lowess <- lowess(df2$alt, f = 0.02)$y
  df2$hr_lowess <- lowess(df2$hr, f = 0.02)$y

  df2$dds0_m  <- c(0, diff( df2$ds0_m ))

  df2$ds0_km  <- df2$ds0_m / 1000
  df2$v0_kmh  <- c(0, diff(df2$ds0_m) ) / df2$dt  * 3.6
  df2$v0_lowess  <-   av_speed(df2$ds0_km, df2$cumtime, nsmooth) #km, h => km/h 


  return( df2 )
}

.4. Загрузка данных. Первые подготовительные операции закончены, теперь можно загрузить реальные данные. Их удобно поместить в какую-то папку, путь к которой будет в переменной data_folder. Я использую для иллюстрации свои поездки по трассам Chess Park.

data_folder <- "C:\\Users\\Vlad\\Documents\\R\\GPS\\data\\"

df1 <- load_tcx( paste0( data_folder, "20190602_094835_Butovo_ChessPark.tcx" ))

nopt <- opt_nsmooth(df1)
nopt

## [1] 26

df2 <- calc_spd(df1, nopt)

.5. Возникает вопрос учитывать ли подъем по высоте для анализа скорости и пройденного пути? Сетевые сервисы типа Endomondo и Strava пишут, что для оценок пройденного пути и скорости они считают Землю плоской, т.е. полностью игнорируют высоту.

Например Strava пишет: “Since this is a GPS-calculated distance, a flat surface is assumed, and vertical speed from topography is not accounted for.
Cons: A flat surface is assumed, and vertical speed from topography is not accounted for. Similar to the above, straight lines connect the GPS points.”

Второй момент: разные методы оценки дистанции — сервисом Endomondo и приближением по эллипсоиду — дают заметно (на 12%) отличающиеся результаты. Более того, если загрузить трек GPX в Strava (этот формат не содержит “подсказок” о пройденном пути, которые есть в формате TCX), то дистанция получится меньше, чем в Endomondo, для одного из моих треков Strava насчитала 14.8км, а Endomondo — 16.6км.

n<- dim(df2)[2]
print( paste0("Endomondo dist./calculated dist.:", df2[n, "ds0_m"]/df2[n,"distance"]/1000))

## [1] "Endomondo dist./calculated dist.:1.1228486352278"

Это расхождение, похоже, систематическое, возможно Endomondo использует для расчета пройденного пути какие-то датчики смартфона, что в целом особенно чувствительно для трасс кросс-кантри, в которых за пару секунд могут случаться заметные движения по вертикали при небольшом движении в горизонтальной плоскости. Трудно сказать. Если мы построим график, на котором длины GPS-отрезков, рассчитанные по геоэллипсоиду, будут отложены против длин, расчитанных Endomondo, то увидим, что есть две группы отрезков — для одних длины по обоим методам приблизительно совпадают, а для других Endomondo рапортует большую длину, чем видно по координатам GPS.

ggplot(df2, aes(x=dds0_m, y = ds))+geom_point()+
  ylab("Geoellipsoid dist., m")+xlab("Reported by Endomondo,m")+ggtitle("Endomondo vs. geoellipse calculated distances")

Эти точки разбросаны равномерно, потому зависимость между этими расстояниями (а потому и скоростями) линейна.

ggplot(df2, aes(x=ds0_m/1000, y=distance_alt))+
  geom_point()+
  ggtitle("Distance, measured by Endomondo vs. Calculated")+
  xlab("Measured by Endomondo, km")+
  ylab("Calculated, km")

Раз уж я загрузил треки в R, то нет смысла повторять то, что и так делают другие. Поэтому дальше я буду пользоваться только расчитанными своими силами скоростями, и с учетом подъема (хотя в среднем разница с “плоской” скоростью получается незначительная).

.6. Для разминки посмотрим график подъемов и спусков. Построение графиков я буду инкапсулировать в специальные функции, которые можно будет потом использовать для других данных.

На рисунке видны как и длинные подъемы, так и быстрые изменения положения по высоте, которые частично вызваны GPS-шумом, а частично соответствуют ныркам на кросс-кантри трейле.

# Plot elevations and smoother

plot_elevations <- function ( df_in ) {
  plot(df_in$cumtime*60, df_in$alt, type = "l", bty = "n",  ylab = "Elevation", xlab = "Time, min", col = "grey40")
  lines(df_in$cumtime*60, df_in$alt_lowess, col = "red", lwd = 3)
  legend(x="bottomright", legend = c("GPS elevation", "LOWESS elevation"),
         col = c("grey40", "red"), lwd = c(1,3), bty = "n")
 title("Elevation vs. time")
}

plot_elevations(df2)

.7. Затем — графики скорости. Их можно строить как относительно времени, так и пройденного расстояния. Как видно из графиков ниже, значительное время своей тренировки я отдыхал, так что средняя скорость за всю тренировку заметно меньше средней скорости на участках, на которых я двигался.

plot_speeds_vs_time <- function( df_in ) {
  n <- dim(df_in)[2]
  plot(df_in$cumtime*60, df_in$speed_kmh_alt, type = "l", bty = "n",  ylab = "Speed (km/h)", xlab = "Time, min",     col = "grey40", main="Speed vs. Time ")
  lines(df_in$cumtime*60, df_in$speed_lowess_alt, col = "blue", lwd = 3)
  legend(x="topright", legend = c("GPS speed", "Smoothed speed"),
         col = c("grey40", "blue"), lwd = c(1,3), bty = "n")
  abline(h = df_in[n, "distance"]/df_in[n,"cumtime"], lty = 2, col = "red")
  abline(h = df_in[n, "distance_alt"]/df_in[n,"cumtime"], lty = 2, col = "yellow")
}
plot_speeds_vs_time( df2 )

# Plot speeds and smoother
plot_speeds_vs_dist <- function( df_in ) {
  n <- dim(df_in)[2]
  plot(df_in$distance, df_in$speed_kmh_alt, type = "l", bty = "n",  ylab = "Speed (km/h)", xlab = "Distance, km",     col = "grey40", main="Speed vs. Distance ")
  lines(df_in$distance, df_in$speed_lowess_alt, col = "blue", lwd = 3)
  legend(x="topright", legend = c("GPS speed", "Smoothed speed"),
         col = c("grey40", "blue"), lwd = c(1,3), bty = "n")
  abline(h = df_in[n, "distance"]/df_in[n,"cumtime"], lty = 2, col = "red")
  abline(h = df_in[n, "distance_alt"]/df_in[n,"cumtime"], lty = 2, col = "yellow")
}
plot_speeds_vs_dist( df2 )

Можно построить гистограмму скоростей, хотя смысла в ней немного.

#Speed histogram
hist_speeds <- function ( df_in ) {
  hist( df_in$speed_kmh, breaks=50, freq= FALSE, xlab="Speed, km/h", main="Distribution of speeds")
}
hist_speeds( df2 )

.8. Теперь перейдем к построению треков. Вначале можно построить просто без всякой карты местности. Для разнообразия цвет каждой точки соответствует скорости передвижения, чем темнее — тем медленнее, чем краснее — тем быстрее.

# Plot the track without any map, the shape of the track is already visible.
plot_track <- function( df_in, maxspeed = 50 ) {
    x1 <- df_in$speed_lowess / maxspeed
  c1 <- rgb( x1, 0, 1-x1 )

  plot(rev(df_in$lon), rev(df_in$lat),  col = c1, lwd = 3, bty = "n", 
       ylab = "Latitude", xlab = "Longitude", pch=20)
}
plot_track( df2 )

.9. Затем — более элегантное построение при помощи Google Maps. Для этого вначале мы вычисляем гео-границы нашего трека, потом получаем для этих границ карту из Google, и при помощи функции ggmap отрисовываем ее. Масштаб карты регулируется параметром zoom. Максимальной детализации соответствует zoom=18. Цвет точек трека будет соответсвовать средней скорости за 90 секунд.

plot_track_on_gmaps <- function( df_in  ) {  

  ###Draw track with google maps

  # mapImageData1 <- get_map(location = c(lon = mean(df1$lon), lat = mean(df1$lon)),
  #                          color = "color",
  #                          source = "google",
  #                          maptype = "satellite",
  #                          zoom = 11)


  sscale = 1
  pad1 <- 2e-4

  box1 <- c( left= min(df_in$lon, na.rm = TRUE) -pad1,
             bottom = min(df_in$lat, na.rm = TRUE)-pad1,
             right = max(df_in$lon,  na.rm = TRUE) +pad1,
             top = max(df_in$lat,  na.rm = TRUE) + pad1)


  mapImageData1 <- get_map(location = box1,
                           color = "color",
                           source = "google",
                           maptype = "satellite", zoom=12 )


  plot1 <- ggmap(mapImageData1,   extent = "panel", #"device"
                 legend="right", padding=0.1, maprange = TRUE) + 
    labs(y = "Latitude", x = "Longitude") + 
    geom_path(aes(x = lon, y = lat, color=speed_lowess), data = df_in,
               size = 2,  show.legend = T) + 
    theme(legend.position = "right")  + 
    scale_color_gradient(low="blue", high="red") 
  print( plot1 )
}
plot_track_on_gmaps( df2)

.10. Так как я регулярно сбиваюсь с пути, то мне любопытно было сопоставить маршрут, по которому я ехал с тем, по которому планировал ехать. Расхождения, прямо скажем, заметные. Планировал ехать по простеньким зелененьким трассам, а по факту форсировал самые сложные в этом парке.

df_plan <- load_gpx( paste0( data_folder, "138pop_2019_by_rockfox.gpx" ))

plot_fact_vs_plan_on_gmaps <- function( df_in, df_plan ) {  

###Draw track with google maps

# mapImageData1 <- get_map(location = c(lon = mean(df1$lon), lat = mean(df1$lon)),
#                          color = "color",
#                          source = "google",
#                          maptype = "satellite",
#                          zoom = 11)


  sscale = 1
  pad1 <- 2e-4

  box1 <- c( left= min(df_in$lon, df_plan$lon, na.rm = TRUE) -pad1,
             bottom = min(df_in$lat, df_plan$lat, na.rm = TRUE)-pad1,
             right = max(df_in$lon, df_plan$lon, na.rm = TRUE) +pad1,
             top = max(df_in$lat,  df_plan$lat, na.rm = TRUE) + pad1)


  mapImageData1 <- get_map(location = box1,
                           color = "color",
                           source = "google",
                           maptype = "satellite", zoom=12 )


  plot1 <- ggmap(mapImageData1,   extent = "panel", #"device"
                 legend="right", padding=0.1, maprange = TRUE) + 
    labs(y = "Latitude", x = "Longitude") + 
    geom_path(aes(x = lon, y = lat, color=speed_lowess), data = df_in,
               size = 1, pch = 20,  show.legend = T) + 
    geom_path(aes(x = lon, y = lat), data = df_plan,
               size = 1, pch = 20,  show.legend = T, , color="yellow") + 
    theme(legend.position = "right")  + 
    scale_color_gradient(low="blue", high="red") 
  print( plot1 )
}
plot_fact_vs_plan_on_gmaps(df2, df_plan)

.11. Интересно глянуть, какой был пульс на тренировке, но в этот раз я ездил без пульсомера, хотя так бывает не всегда. Потому мне потребуется другой трек. Принципы построения графиков все те же, так что я на них не буду останавливаться.

df3 <-  load_tcx( paste0( data_folder, "20190529_115118.tcx" ))
nopt <- opt_nsmooth(df3)
nopt

## [1] 48

df3 <- calc_spd(df3)

plot_hr_vs_time <- function( df_in ) {
  plot(df_in$cumtime*60, df_in$hr, type = "l", bty = "n",  ylab = "HR (bpm)", xlab = "Time, min",  
       col = "grey40", main="HR vs. Time ")
  lines(df_in$cumtime*60, df_in$hr_lowess, col = "blue", lwd = 3)
  legend(x="topright", legend = c("HR", "Averaged HR"),
         col = c("grey40", "blue"), lwd = c(1,3), bty = "n")
  abline(h = mean(df_in$hr), lty = 2, col = "red")
  abline(h = mean(df_in$hr_lowess), lty = 2, col = "yellow")
}
plot_hr_vs_time( df3 )

plot_hr_vs_dist <- function( df_in ) {
  plot(df_in$distance, df_in$hr, type = "l", bty = "n",  ylab = "HR (bpm)", xlab = "Distance, km",  
       col = "grey40", main="HR vs. Distance ")
  lines(df_in$distance, df_in$hr_lowess, col = "blue", lwd = 3)
  legend(x="topright", legend = c("HR", "Averaged HR"),
         col = c("grey40", "blue"), lwd = c(1,3), bty = "n")
  abline(h = mean(df_in$hr), lty = 2, col = "red")
  abline(h = mean(df_in$hr_lowess), lty = 2, col = "yellow")
}
plot_hr_vs_dist( df3 )

plot_hr_on_ggmap <- function( df_in ) {
  sscale = 1
  pad1 <- 2e-4

  box1 <- c( left= min(df_in$lon, na.rm = TRUE) -pad1,
             bottom = min(df_in$lat, na.rm = TRUE)-pad1,
             right = max(df_in$lon, na.rm = TRUE) +pad1,
             top = max(df_in$lat,  na.rm = TRUE) + pad1)


  mapImageData1 <- get_map(location = box1,
                           color = "color",
                           source = "google",
                           maptype = "satellite", zoom=12 )


    ggmap(mapImageData1,   extent = "panel", #"device"
                 legend="right", padding=0.1, maprange = TRUE) + 
    labs(y = "Latitude", x = "Longitude") + 
    geom_path(aes(x = lon, y = lat, color=hr_lowess), data = df_in,
               size = 2, show.legend = T) + 
    theme(legend.position = "right")  + 
    scale_color_gradient(low="blue", high="red") +
    ggtitle("Pulse rate along the track")
}
plot_hr_on_ggmap(df3)

.12. Ну что же, пульс у меня высокий. И чем сложнее участок трассы, тем он выше. На обычных покатушках по лесу, без выезда на трейлы красного уровня, он в среднем около 130, с редкими заходами на 150. Посмотрим теперь, сколько времени в каких пульсовых зонах я провожу. Ну в общем…нехорошо.

hr_zone  <- function(hr, hrmax) {
  return( as.integer( floor( ( hr/hrmax - 0.5 ) * 10 ) + 1 ))
}

plot_hr_stat <- function( df_in, hrmax = 175 ) {
  x1 <-  df_in[  , .(dtsum = sum(dt)), by =.(hr)]
  setkey(x1, hr)

  gg1<- ggplot(x1, aes(x=hr, y=dtsum/60))+
    geom_line(color="red")+
    ggtitle("Time spent at given HR")+
    geom_smooth() +
    xlab("Heart rate, bmp")+ 
    ylab("Time spent, min") 

  gg_add_zones <- function(gg1, hrmax) { 
    for( i in c(2:6)) {
      h <- hrmax * (0.4 + i* 0.1)
      gg1 <- gg1 + geom_vline(xintercept =  h, linetype= "dotted")
    }

    i <- c(2:5)
    df_labels <- data.frame( hr = hrmax * (0.45 + i* 0.1), hr_zone = i )

    gg1 <- gg1 + geom_text( data=df_labels, aes( x=hr, 0, label=paste0( "Zone ", hr_zone), hjust = 0.5 ))
  }
  gg1 <- gg_add_zones( gg1, hrmax )

  print(gg1)

  gg2<- ggplot(x1, aes(x=hr, y=dtsum/sum(dtsum)))+
    geom_line(color="red")+
    geom_smooth() +
    ggtitle("Share of time spent at given HR")+
    xlab("Heart rate, bmp")+ 
    ylab("Share of time")

  gg2 <- gg_add_zones( gg2, hrmax )
  print(gg2)

  gg3 <- ggplot(x1, aes(x=hr, y=cumsum(dtsum)/sum(dtsum)))+
    geom_line(color="red") +
    ggtitle("Share of time spent at HR<x")+
    xlab("Heart rate, bmp")+ ylab("Share of time")
  gg3 <- gg_add_zones( gg3, hrmax )

  print(gg3)

  gg4 <- ggplot(x1, aes(x=hr, y=cumsum(dtsum)/60))+
    geom_line(color="red") +
    ggtitle("Time spent at HR<x")+
    xlab("Heart rate, bmp")+ ylab("Time,min")
  gg4 <- gg_add_zones( gg4, hrmax )

  print(gg4)

  }

plot_hr_stat(df3)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

.13. Пульсовые зоны можно также отрисовать на треке и сопоставить с участками трассы. Ожидаемо, чем сложнее для меня участок, тем выше пульс.

df3$hr_zone <- as.integer( hr_zone(df3$hr, hrmax))  

plot_hrzone_on_ggmap <- function( df_in ) {
  sscale = 1
  pad1 <- 2e-4

  box1 <- c( left= min(df_in$lon, na.rm = TRUE) -pad1,
             bottom = min(df_in$lat, na.rm = TRUE)-pad1,
             right = max(df_in$lon, na.rm = TRUE) +pad1,
             top = max(df_in$lat,  na.rm = TRUE) + pad1)


  mapImageData1 <- get_map(location = box1,
                           color = "color",
                           source = "google",
                           maptype = "satellite", zoom=12 )


  ggmap(mapImageData1,   extent = "panel", #"device"
        legend="right", padding=0.1, maprange = TRUE) + 
    labs(y = "Latitude", x = "Longitude") + 
    geom_point(aes(x = lon, y = lat, color=as.factor(hr_zone)), data = df_in,
              size = 2, show.legend = T) + 
    theme(legend.position = "right")  + 
    ggtitle("Pulse zones on the track") +
    scale_color_discrete()
}
plot_hrzone_on_ggmap(df3)

“Шестая зона” особого смысла не имеет и означает, что часть времени, хоть и совсем ничтожную, я провожу с пульсом, выше расчетно-предельного по возрасту.