From a GPX file to research-grade trail measures

This vignette walks through the metrictrails pipeline with the files shipped in the package: a trail run on Manquehue hill (Santiago, Chile) recorded by a GNSS watch, the same loop drawn as a route, and a small Digital Elevation Model (DEM) of the area.

library(metrictrails)

activity_file = system.file("extdata", "gpx_activity.gpx", package = "metrictrails")
route_file = system.file("extdata", "gpx_route.gpx", package = "metrictrails")
dem_file = system.file("extdata", "test_dem.tif", package = "metrictrails")

The pipeline has five steps:

  1. ingest a GPX file, read_stream();
  2. apply the naming convention, uniform_stream();
  3. sanitize the coordinates, sanitize_stream();
  4. resample to a uniform spatial resolution, resample_stream();
  5. replace the elevation with values from a DEM, ele_correction().

1. Ingest

read_stream() reads the track points of a GPX file with sf. Every track segment becomes a value of segment, and Garmin point extensions (heart rate, cadence, temperature) become numeric columns.

raw = read_stream(activity_file)
raw
#> # A tibble: 2,510 × 8
#>      lon   lat   ele timestamp           segment atemp    hr   cad
#>    <dbl> <dbl> <dbl> <dttm>                <int> <dbl> <dbl> <dbl>
#>  1 -70.6 -33.4  772  2020-08-10 19:36:38       1    27   106    47
#>  2 -70.6 -33.4  771  2020-08-10 19:36:39       1    27   108    47
#>  3 -70.6 -33.4  772. 2020-08-10 19:36:42       1    27   111    76
#>  4 -70.6 -33.4  771. 2020-08-10 19:36:43       1    27   114    76
#>  5 -70.6 -33.4  771. 2020-08-10 19:36:44       1    27   117    77
#>  6 -70.6 -33.4  771. 2020-08-10 19:36:47       1    27   121    79
#>  7 -70.6 -33.4  771. 2020-08-10 19:36:49       1    27   125    79
#>  8 -70.6 -33.4  771. 2020-08-10 19:36:50       1    27   127    79
#>  9 -70.6 -33.4  772. 2020-08-10 19:36:53       1    27   130    79
#> 10 -70.6 -33.4  772. 2020-08-10 19:36:56       1    27   133    80
#> # ℹ 2,500 more rows

2. Naming convention

uniform_stream() renames the columns to lon, lat, ele, timestamp and time (seconds since the start). It also works on data frames from other sources, such as the Strava API streams in strava_streams.

raw = uniform_stream(raw)
names(raw)
#> [1] "lon"       "lat"       "ele"       "timestamp" "time"      "segment"  
#> [7] "atemp"     "hr"        "cad"
head(uniform_stream(strava_streams), 3)
#> # A tibble: 3 × 13
#>     lon   lat   ele  time cadence distance grade_smooth heartrate moving  temp
#>   <dbl> <dbl> <dbl> <dbl>   <int>    <dbl>        <dbl>     <int> <lgl>  <int>
#> 1 -70.6 -33.4 1052.     0      NA   0              10.7        NA FALSE     NA
#> 2 -70.6 -33.4 1053.    10      NA   0.0079          7.5        NA TRUE      NA
#> 3 -70.6 -33.4 1054.    20      NA   0.0112          3.1        NA FALSE     NA
#> # ℹ 3 more variables: velocity_smooth <dbl>, id <dbl>, type <chr>

3. Sanitize

sanitize_stream() first flags readings that are physically implausible: speed, acceleration, jumps at the edges, slope, vertical speed, and elevation changes while the athlete stands still. Then it smooths lon and lat. The default method is a constant-velocity Kalman filter with a Rauch-Tung-Striebel smoother; LOESS, a Gaussian kernel, Savitzky-Golay and the full post-processing pipeline of Martinez, Sanchez & Diaz (2026) are also available.

clean = sanitize_stream(raw)
attr(clean, "sanitize_report")$counts
#>       flag_time      flag_speed      flag_accel       flag_jump      flag_slope 
#>               0               0               0               0              71 
#> flag_vert_speed   flag_idle_ele        flag_gap    flag_low_res 
#>               3               3               0              10

The flags are computed on the raw readings, before smoothing. In this activity the positions raise no flags. The flags concern the raw elevation: slopes above 50 degrees between close readings, implausible vertical speeds, and elevation changes while standing still. Flagged readings can be kept (action = "flag", the default), dropped or interpolated.

methods = c("kalman", "loess", "gaussian", "savgol", "pipeline")
lengths = sapply(methods, function(m) {
  summarise_stream(sanitize_stream(raw, method = m))$distance_2d
})
round(c(raw = summarise_stream(raw)$distance_2d, lengths))
#>      raw   kalman    loess gaussian   savgol pipeline 
#>     7847     7452     6253     6686     7454     7523

Every method removes part of the positional noise, which shortens the path. LOESS with its default span smooths the most.

4. Resample: the Nyquist resolution

A watch that records every second produces a spatial resolution that is a random variable: the distance between readings is the speed times the sampling interval. Because trails are fractal, distance and elevation gain depend on that resolution (Sanchez et al., 2025), so two recordings of the same trail are not directly comparable.

resolution_stream() describes the distribution of the horizontal distance between moving readings:

resolution_stream(clean)
#> # A tibble: 1 × 12
#>       n pct_idle  mean    sd   q05   q25   q50   q75   q95   max     q   res
#>   <int>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1  2510   0.0797  2.97  3.07 0.512  1.08  2.00  3.49  9.10  22.8   0.9  6.31

Its 0.9 quantile, res, is the finest uniform spacing that represents the path without inventing readings between sparse fixes. resample_stream() uses it by default:

resampled = resample_stream(clean)
attr(resampled, "res")
#> [1] 6.306554
nrow(resampled)
#> [1] 1182

The fractal effect is easy to see by resampling the same stream at increasing resolutions:

res = c(1, 2, 3, 5, 8, 13, 20, 30, 50, 80)
distance = sapply(res, function(r) summarise_stream(resample_stream(clean, res = r))$distance_2d)
plot(res, distance / 1000, log = "x", type = "b", pch = 19,
     xlab = "Resolution (m, log scale)", ylab = "Distance (km)")
abline(v = attr(resampled, "res"), lty = 2)

To compare several streams, resample all of them with one common resolution, at least the largest of their Nyquist resolutions:

streams = list(activity = clean, route = read_stream(route_file))
resolution_stream(streams)[, c("stream", "q50", "res")]
#> # A tibble: 2 × 3
#>   stream     q50   res
#>   <chr>    <dbl> <dbl>
#> 1 activity  2.00  6.31
#> 2 route    13.0  37.2
common = max(resolution_stream(streams)$res)

The route was drawn with sparse vertices, so it sets the common resolution.

fractal_stream() summarises the same behaviour as a fractal dimension (divider method):

fd = fractal_stream(clean)
fd$dimension
#> [1] 1.073242

5. Elevation from a DEM

Raw elevation comes from a barometer or from GNSS, and it is noisy and often rounded to the metre. ele_correction() replaces it with values from a DEM. Without a dem argument, get_dem() downloads AWS Terrain Tiles at zoom 14 along a corridor around the stream. Here we use the DEM shipped with the package:

final = ele_correction(resampled, dem = dem_file)
# same result with a downloaded DEM, refined to 1 m by bilinear interpolation
dem = get_dem(resampled, res = 1)
final = ele_correction(resampled, dem = dem)
d_raw = differential_stream(raw)
d_final = differential_stream(final)
plot(d_raw$distance / 1000, raw$ele, type = "l", col = "grey60",
     xlab = "Distance (km)", ylab = "Elevation (m)")
lines(d_final$distance / 1000, final$ele, lwd = 1.5)
legend("topright", c("device", "pipeline + DEM"), col = c("grey60", "black"), lwd = c(1, 1.5), bty = "n")

Measures

differential_stream() adds per-reading differences and rates, and summarise_stream() the course-level measures, here before and after the pipeline:

cols = c("distance_2d", "elevation_gain", "elevation_loss", "km_effort", "resolution")
cbind(stream = c("raw", "pipeline"),
      rbind(summarise_stream(raw)[, cols], summarise_stream(final)[, cols]))
#>     stream distance_2d elevation_gain elevation_loss km_effort resolution
#> 1      raw    7847.284       889.2000       875.4000  16.73928   3.127654
#> 2 pipeline    7280.702       950.9853       945.9306  16.79056   6.164862

Tobler parameters

tobler_stream() fits an asymmetric Tobler hiking function to the speed-slope relation. It models the frontier of performance (by default the 0.9 quantile of speed), not the mean, by minimising the quantile loss.

params = tobler_stream(final)
params
#> # A tibble: 1 × 8
#>      v1    m1 c_plus c_minus   tau     n    r1 method 
#>   <dbl> <dbl>  <dbl>   <dbl> <dbl> <int> <dbl> <chr>  
#> 1  3.30 0.148   3.40 -0.0248   0.9  1083 0.236 pinball

d = differential_stream(final)
keep = is.finite(d$speed) & is.finite(d$slope) & abs(d$slope) <= 0.6
plot(d$slope[keep], d$speed[keep], pch = 16, cex = 0.4, col = "grey60",
     xlab = "Slope (rise over run)", ylab = "Horizontal speed (m/s)")
slope = seq(-0.6, 0.6, by = 0.01)
lines(slope, tobler_speed(slope, params), lwd = 2)

In this single activity the frontier is flat downhill (c_minus close to 0) and the speed only drops on steep climbs. Many activities, or a population of athletes, give a better-defined frontier.

Two streams

detect_stream() is a predicate: does the activity follow the route, in which direction, and between which readings?

route = read_stream(route_file)
detect_stream(raw, route)
#> # A tibble: 1 × 6
#>    pass start_row end_row direction coverage status
#>   <int>     <int>   <int> <chr>        <dbl> <chr> 
#> 1     1         1    2506 forward          1 match

locate_stream() projects every reading on the route (linear referencing), and compare_stream() measures the spatial agreement of two streams resampled at the same resolution:

located = locate_stream(clean, route)
summary(abs(located$offset_cross))
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#>  0.01153  3.30181  5.52755  6.92778  8.84520 39.74600

compare_stream(resample_stream(clean, res = 2), resample_stream(route, res = 2))
#> # A tibble: 1 × 8
#>    dims hausdorff hausdorff_ab hausdorff_ba separation_ab separation_ba
#>   <dbl>     <dbl>        <dbl>        <dbl>         <dbl>         <dbl>
#> 1     2      40.0         40.0         33.9          5.28          5.10
#> # ℹ 2 more variables: metric_score_ab <dbl>, metric_score_ba <dbl>

Matching to a known reference

When the athlete is known to have followed a reference trail, for example rounds on a trail measured with a high-precision protocol, match_stream() estimates where the athlete was along the reference at every timestamp. The match never moves backwards and never exceeds the Tobler ceiling, and the global optimum is found by dynamic programming. The result keeps the time and the other variables of the stream, takes the coordinates and elevation of the reference, and reports the offset of every reading along and across the trail.

reference = route %>% resample_stream(res = 2) %>% ele_correction(dem = dem_file)
matched = match_stream(clean[1:800, ], reference, params = params)
attr(matched, "match_report")
#> # A tibble: 1 × 10
#>   group direction     n vertices_used max_per_vertex vertices_over_cap
#>   <lgl> <chr>     <int>         <int>          <int>             <int>
#> 1 NA    forward     800           497              8                72
#> # ℹ 4 more variables: pct_at_ceiling <dbl>, objective <dbl>,
#> #   offset_median <dbl>, offset_q95 <dbl>
summary(matched[, c("offset_d", "offset_along", "offset_cross")])
#>     offset_d        offset_along       offset_cross    
#>  Min.   : 0.2904   Min.   :-23.9252   Min.   :-30.416  
#>  1st Qu.: 4.7058   1st Qu.: -1.0220   1st Qu.: -8.146  
#>  Median : 7.4410   Median :  0.2385   Median : -5.060  
#>  Mean   : 8.4524   Mean   :  0.8401   Mean   : -5.630  
#>  3rd Qu.:11.2321   3rd Qu.:  2.9667   3rd Qu.: -2.649  
#>  Max.   :30.8838   Max.   : 24.5885   Max.   : 10.217
plot(matched$time / 60, matched$s_ref / 1000, type = "l",
     xlab = "Time (min)", ylab = "Position along the reference (km)")

Writing the result

out = tempfile(fileext = ".gpx")
write_stream(final, out, name = "Manquehue loop")
nrow(read_stream(out))
#> [1] 1182

References