Source files of https://rpubs.com/aoles/cyclingWeightings

coordinates_bayern-cycling-1000.Rda

cyclingWeightings.Rmd

---
title: "Cycling weightings compared"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Preprocessing of data

```{r}
library(openrouteservice)
library(leaflet)

load("recommendended-bike-1000.Rda")

n = 1000

lapply(res, function(routes) {
  sapply(routes, function(x) x$features[[1]]$properties$summary$distance/1000) 
}) %>% data.frame -> distance

lapply(res, function(routes) {
  sapply(routes, function(x) x$features[[1]]$properties$summary$duration)
}) %>% data.frame -> duration

lapply(res, function(routes) {
  sapply(routes, function(x) x$features[[1]]$properties$ascent) 
}) %>% data.frame -> ascent

lapply(res, function(routes) {
  lapply(routes, function(x) {
    x = x$features[[1]]$properties$extras$suitability$summary
    values = sapply(x, `[[`, "value") 
    amount = sapply(x, `[[`, "amount")
    v = rep(0, 8)
    v[values-2] <- amount
    v
  })
}) %>% lapply(function(routes) {
  sapply(routes, function(v) {
    sum(0:7 * v/100)
  })
}) %>% data.frame -> suitability
```


## Route characteristics

```{r}
range(distance$fastest)
```

```{r}
hist(distance$fastest, xlim = c(0, 500), breaks = 50)
```


## Suitability

```{r}
compare <- function(what, ref, f) {
  sapply(lapply(what[setdiff(names(what), ref)], f, what[ref]), sum)
}
```

Number of routes that have a better suitability index compared to GraphHopper's `fastest`...

```{r}
compare(suitability, "fastest", `>`)
```

.. and a worse suitability.

```{r}
compare(suitability, "fastest", `<`)
```



Even though current `recommended` is a better choice than `fastest`, all the other alternatives seem to perform better with the proposed `new` being the best.

But how do actually compare the other weightings to `new`? 
```{r}
compare(suitability, "new", `<`)
compare(suitability, "new", `>`)
```

Routes which become worse with the new weighting.

```{r}
(worse_suitability <- which(suitability$new < suitability$fastest | suitability$new < suitability$recommended))
```

```{r, out.width = '100%', out.height = '780px', echo = FALSE}
leaflet() %>% 
  addTiles(group = "OpenStreetMap") %>%
  addTiles("https://dev.{s}.tile.openstreetmap.fr/cyclosm/{z}/{x}/{y}.png", group ="CyclOSM" ) %>%
  addTiles("https://{s}.tile.thunderforest.com/cycle/{z}/{x}/{y}.png?apikey=13efc496ac0b486ea05691c820824f5f", group ="OpenCycleMap") %>%
  addGeoJSON(res$fastest[worse_suitability], fill=FALSE, color = c("#000"), group = "current fastest") %>%
  addGeoJSON(res$recommended[worse_suitability], fill=FALSE, color = c("#f00"), group = "current recommended") %>%
  addGeoJSON(res$new[worse_suitability], fill=FALSE, color = c("#0f0"), group = "new") %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap", "CyclOSM", "OpenCycleMap"),
    overlayGroups = c("current fastest", "current recommended", "new"),
    position = "bottomleft"
  ) %>%
  fitBBox(c(8.974646, 47.265430, 13.849470, 50.567790))
```


## Ascent

Total ascent values generated across routes.

```{r}
sapply(ascent, sum)
```

`new` is second best even though it produces the longest routes.

```{r}
sapply(distance, sum)
```

Number of routes that have a lower ascent (better) compared to GraphHopper's `fastest` ...

```{r}
compare(ascent, "fastest", `<`)
```

... and higher ascent (worse).

```{r}
compare(ascent, "fastest", `>`)
```

Again, the `new` weighting seems to be the best alternative to `fastest`.

## Summary

Percentage of `recommended` routes which are worse than `fastest` in terms of both suitability and ascent.

```{r}
sum(ascent$recommended > ascent$fastest & suitability$recommended < suitability$fastest) / n * 100
```

Percentage of `new` routes which are worse than `fastest` in terms of both suitability and ascent.

```{r}
sum(ascent$new > ascent$fastest & suitability$new < suitability$fastest) / n * 100
```

Percentage of `new` routes which are worse than `recommended` in terms of both suitability and ascent.

```{r}
sum(ascent$new > ascent$recommended & suitability$new < suitability$recommended)  / n * 100
```

generate_routes.R

library("openrouteservice")

set.seed(0L)
options(openrouteservice.url = "http://localhost:8082/ors")

profile = "cycling-regular"

n = 1000
rep = 1

### generate n random coordinates and convert them to start/endpoints

coordinates_file = "coordinates_bayern-cycling-1000.Rda"

if (file.exists(coordinates_file)) {
  load(coordinates_file)
} else {
  library("sf")
  library("geojson")
  
  x = readLines("~/Data/polygons/bayern.geojson")
  bbox = geo_bbox(x)
  poly = st_read(paste(x, collapse = ""))
  
  pts = vector(mode = "list", n)
  
  i = 0;
  while (i < n) {
    coords = c(runif(1, bbox[1], bbox[3]), runif(1, bbox[2], bbox[4]))
    if (st_within(st_point(coords), poly, sparse=FALSE))
      if ( !inherits(try(ors_directions(list(coords, coords), profile=profile), silent=TRUE), "try-error") )
        pts[[(i = i+1)]] <- coords
  }
  
  # duration matrix with diagonal excluded
  res <- ors_matrix(pts, profile = profile)
  durations <- res$durations
  durations[seq(1, by = n+1, length.out = n)] <- NA
  max <- max(res$durations)
           
  d <- vector(mode = "numeric", n)
  coordinates <- vector(mode = "list", n)
           
  for (i in 1:n) {
    r <- runif(1, max = max)
    w <- which.min(abs(durations - r))
    d[i] <- durations[w]
    index <- arrayInd(w, dim(durations))
    a = index[1,1]
    b = index[1,2]
    durations[a, b] <- NA
    durations[b, a] <- NA
    coordinates[[i]] <- list(pts[[a]], pts[[b]])
  }
  
  coordinates <- coordinates[order(d)]
  
  save(coordinates, file = coordinates_file, compress = "xz")
}

profile_args = list(profile = profile, format = 'geojson', instructions = FALSE, elevation = TRUE, extra_info = "suitability")

profile_args = list(
  fastest = c(profile_args, preference = "fastest"),
  recommended = c(profile_args, preference = "recommended"),
  reduced = c(profile_args, preference = "recommendedr"),
  scaled = c(profile_args, preference = "recommendeds"),
  new = c(profile_args, preference = "recommendednew")
  )

res = lapply(names(profile_args), function(name) {
  lapply(1:n, function(id) {
    cat("                         \r", name, id)
    label = sprintf("[BENCHMARK] %s; %d", name, id)
    cl = as.call(c(ors_directions, coordinates[id], profile_args[[name]], id = label))
    res = try(eval(cl), silent = TRUE)
    if ( inherits(res, "try-error") ) {
      attr(res, "condition")$message
    } else {
      cl$id = NULL
      query_times = replicate(rep-1, attr(eval(cl), "query_time"))
      attr(res, "query_time") = c(attr(res, "query_time"), query_times)
      res
    }
  })
})

names(res) <- names(profile_args)

save(res, file = "recommendended-bike-1000.Rda", compress = "xz")

recommendended-bike-1000.Rda