Travel time calculation with R and data visualization with Observable

Application to artifical climbing walls in Paris and its neighbourhood

Author
Affiliation
Ronan Ysebaert

UAR RIATE, Université Paris Cité, CNRS

Published

January 20, 2023

The aim of this notebook consists in showing how to build, visualize and reproduce accessibility indicators by combining points of interest (POI) coming from the OpenStreetMap database, socio-economic indicators included in small territorial division (IRIS) coming from institutional data source (INSEE) and routing engines (OSRM).

The graphical outputs displayed in this notebook are further developed in an Observable collection. To introduce the reader to the issues raised by indoor sport-climbing in Paris, have a look to this notebook.

This Quarto document combines 2 programming languages : R for data processing, and ObservableJS for data visualizations.

Data processing with R

The entire R script can be found here. Geojson files resulting from the data processing are available in the data-conso folder of the github repository of the project. These files are directly imported in Observable notebooks for data visualization.

Data sources

Three data sources, coming from three data providers, will be used:

Input data is not provided in the GitHub repository (too large files), but can easily be downloaded and used (no constraints of use).

Map layout

A bounding box of 5 km around Paris, without Bois de Boulogne and Bois de Vincennes for a better centering of the map layout around Paris is created.

A second bounding box (10 km around Paris) is created to catch OpenStreetMap POI and IRIS in the neighbourhood of this study area and avoid “border effects” for the upcoming travel-time indicators calculations.

Code 
# 1. Map layout preparation at IRIS scale (source IGN)----
library(sf)
iris <- st_read("data-raw/CONTOURS-IRIS.shp", quiet = TRUE)

# Extract Paris and delete 
# Bois-de-Vincennes / Boulogne Iris for map template
iris$dep <- substr(iris$INSEE_COM, 1, 2)
paris <- iris[iris$dep == "75",]
paris <- paris[!paris$NOM_IRIS %in% c("Bois de Vincennes 1",
                                      "Bois de Vincennes 2",
                                      "Bois de Boulogne 1",
                                      "Bois de Boulogne 2",
                                      "Bois de Boulogne 3"),]
paris <- st_union(paris)

# 5 km around Paris map layout
paris5k <- st_buffer(paris, 5000)
paris5k <- st_as_sfc(st_bbox(paris5k, crs = 2154))
paris <- paris5k

# 10 km around Paris (get OSM data) in long/lat
paris10k <- st_buffer(paris, 10000)
paris10k <- st_as_sfc(st_bbox(paris10k, crs = 2154))

# Intersection with IRIS
iris10k <- st_intersection(iris, paris10k)

# Bounding box for osm extract
paris10k <- st_transform(paris10k, 4326)
paris10k <- st_bbox(paris10k) 
paris10k <- as.vector(paris10k)

Feed IRIS with socio-economic data (INSEE)

Data is enriched by socio-economic data (disposable median income 2019 and total population 2018) for further analysis. We keep only “Habitation” IRIS (dedicated to dwellings) for origins-destinations calculations.

Code 
# 2. Feed IRIS layer by socio-economic data (INSEE) ----
library(readxl)

# Population 2018 (IRIS)
df <- read_xlsx("data-raw/base-ic-evol-struct-pop-2018.xlsx", skip = 5, sheet = "IRIS")
iris10k <- merge(iris10k[,c("CODE_IRIS", "NOM_IRIS", "TYP_IRIS", "NOM_COM")], 
                 df[,c("IRIS","P18_POP")],
                 by.x = "CODE_IRIS", by.y = "IRIS", all.x = TRUE)

# Median Income (2018)
df <- read_xlsx("data-raw/BASE_TD_FILO_DISP_IRIS_2019.xlsx", skip = 5, sheet = "IRIS_DISP")
iris10k <- merge(iris10k, df[,c("IRIS","DISP_MED19")],
                 by.x = "CODE_IRIS", by.y = "IRIS", all.x = TRUE)

# Intersection with study area
iris <- st_intersection(iris10k, paris5k)

Prepare IRIS for travel-time calculation

IRIS centroids are extracted. These points will be used for the origins of travel-time calculations. Only IRIS including dwellings are kept (TYP_IRIS == “H”).

These sf objects are transformed in latitude/longitude (origin-destination calculations requirements and for final export in geojson format).

Code 
# Keep only habitation IRIS for origins calculation
ori <- iris10k[iris10k$TYP_IRIS == "H",]
ori <- st_centroid(ori)

# Transform in long/lat
ori <- st_transform(ori, crs = 4326)
iris <- st_transform(iris, crs = 4326)

Import OSM points of interest

Points of interest (climbing areas) are extracted from OpenStreetMap thanks to the osmdata R package (Mark Padgham (2022)).

To access to these OSM features, the OSM key-value pair (or OSM tag) must be set. For climbing areas, the most appropriate is sport=climbing. In the OpenStreetMap wiki, it is mentioned that “sport=climbing should be preferably applied to noded for artificial climbing walls”. Thus, we only consider points responding to the query.

The geographical coverage of the query covers a bounding box of 10 km around Paris (paris10k`).

Code 
# 3. Extract OSM objects (climbing and map layout)----
library(osmdata)

# define a bounding box
q0 <- opq(bbox = paris10k) 

# extract climbing areas
q <- add_osm_feature(opq = q0, key = 'sport', value = "climbing")
res <- osmdata_sf(q)
dest <- res$osm_points
dest[,"name"] <- iconv(dest$name, from = "UTF-8", to = "UTF-8")

The OSM attributes are transformed afterwards for further analytical purposes in order to differentiate private and associative structures, and bouldering areas and climbing walls (to different climbing practices).

I argue that the accuracy and completeness of the data is quite good : I have edited missing points on OpenStreetMap database with my personal knowledge and upstream investigation :-)

For the presentation of the artificial climbing landscape in Paris and the difference between private - FSGT and FFME structures, have a look to this notebook.

Code 
# Cleaning
private <- dest[!is.na(dest$brand),] # Manage private and associative areas
asso <- dest[!is.na(dest$federation),]
asso$type <- "Associative structure"
private$type <- "Speculative structure"
dest <- rbind(asso, private)
dest$federation[is.na(dest$federation)] <- "Private"
# Find walls and boulders
dest[c("climbing.toprope", "climbing.boulder")][is.na(dest[c("climbing.toprope", "climbing.boulder")])] <- "no"
dest$climbing_type <- ifelse(dest$climbing.toprope == 'yes' & 
                               dest$climbing.boulder == "yes", 'Wall and bouldering',
                             ifelse(dest$climbing.toprope == 'yes' & 
                                      dest$climbing.boulder == "no" , 'Wall',
                                    ifelse(dest$climbing.toprope == 'no' & 
                                             dest$climbing.boulder == "yes" ,
                                           'Bouldering', NA)))
# Keep only attributes of interest and rename it
cols <- c("osm_id", "name", "climbing_type", "climbing.length",
          "climbing.routes", "type", "federation", "brand")
dest <- dest[,cols]
colnames(dest)[4:5] <- c("climbing_length", "climbing_routes")

# Intersection with bounding box 
poi <- st_transform(dest, 2154)
poi <- st_intersection(poi, paris5k)
poi <- st_transform(poi, 4326)

The data preparation allows to prepare origins-destinations layers for travel-time calculation. Origins correspond to the IRIS layer (only type H : 2141 points). Destinations to artificial climbing areas (45 points).

Code 
library(mapsf)

# Extract centroids from IRIS dedicated to dwellings
ori <- st_centroid(iris)
ori <- ori[ori$TYP_IRIS == "H",]

# Map creation
par(mar = c(0,0,0,0))

# Set a map theme (library mapsf)
my_theme <- list(bg = NA, fg = NA, mar = c(0, 0, 0, 0), tab = TRUE, pos = "left", 
                  inner = TRUE, line = 1.3, cex = 1, font = 2)

# Create the map
mf_init(com, expandBB = c(0,0.4,0,0), theme = my_theme)

mf_typo(
  x = iris, var = "TYP_IRIS", 
  pal = c("peachpuff", "#feb8ff", "#f0e6f0", "#f1f7ab"), lwd = 0.25,
  val_order = c("H", "A", "D", "Z"), border = "white", leg_pos = c(631000, 6872000),
  leg_title = "IRIS types (H = dwellings)", add = TRUE)

mf_map(ori, pch = 20, col = "red", cex = .2, add = TRUE)

mf_legend(
  type = "symb", pos =  c(631000, 6864000), val = c("H", "H"),
  pt_pch = c(20,20), pt_cex =.2, title = "Origins : IRIS centroids type H",
  pal = c("red", "red"))

mf_typo(x = poi, var = "federation", pch = 21, cex = .7,
        val_order = c("Private", "FSGT", "FFME"),
        pal = c("#377eb8", "#e41a1c", "#ff7f00"),
        leg_title = "Destinations:\nArtificial climbing areas",
        leg_pos = c(631000, 6860000),
        add = TRUE)

mf_map(com, col = NA, border = "black", add = TRUE)

mf_title("Layers presentation for origins-destinations calculation")
mf_scale(size = 5, col = "black", pos = "bottomright")
mf_credits(paste0("Sources : © OpenStreetMap and Contributors, IGN, INSEE, 2022\n",
                  "Realisation : Ronan Ysebaert, 2022"))

Travel-time calculations with OSRM

Origins and destinations required to compute travel-time calculations are now available. The computation is realized using the osrm R package (Giraud (2022b)). This package allows the computation of routes, trips, isochrones and travel distances matrices (travel time and kilometer distance).

Considering that the input data are quite important (more than 2000 origin points and 48 destination points) and to avoid to overload the OSRM demo server, I have run my own instance of OSRM based on a docker container solution. The procedure to implement this solution is explained in the following documentation, or more specifically for Windows here.

Once it is done, the connection to the routing engine is operational locally (with the URL http://localhost:5000/` for me).

With OSRM, it is possible to choose several profiles depending on the routing we want to use (car, bike, walking). In our case, we choose the bike profile: Cycling is a widespread practice (with public transport) for climbers to reach their activities in this kind of urban context. Moreover, it allows also to avoid to use the car profile, which underestimates the travel-time in metropolitan areas (traffic congestion not taken into account).

The bicycle profile is described in the github OSRM repository. Basically, the default speed is 15 km/h. It avoids the access to highways, reduce the driving speed by 30 % for unsafe roads. It is important to keep in mind that landforms (elevation) are not considered in the default profile. Some thoughts and solutions exist on the subject, using elevation rasters (Liedman (2022) ; Mapbox (2015)). This will not be considered here.

IRIS to climbing areas

The code below compute travel-time calculation between the geometric centre of IRIS (ori) and artificial climbing areas (dest). It is done for all climbing structures and reproduced according to the type of climbing structure : private or associative structures. For associative structures, we distinguish the FSGT and FFME federations, that do not have exactly the same goals in term of practice.

Resulting travel-time matrixes are exported in the data-conso folder for other possible uses (and to avoid running these important calculations systematically).

Code 
# 4. Origin - Destination calculation with OSRM ----
library(osrm)

# Manage ids
row.names(ori) <- ori$CODE_IRIS

# Connexion to osrm local server (the osrm container must running under docker)
options(osrm.server = "http://localhost:5000/", osrm.profile = "bike")

# Origin-destination calculation
# All structures
df <- osrmTable(src = ori, dst = dest, measure = "duration")
df <- data.frame(df$duration)
colnames(df) <- as.character(dest$osm_id)
row.names(df) <- as.character(ori$CODE_IRIS)
write.csv(df, "data-conso/bike-duration.csv")

# Private structures
dest_priv <- dest[dest$type == "Speculative structure",] 
df2 <- osrmTable(src = ori, dst = dest_priv, measure = "duration")
df2 <- data.frame(df2$duration)
colnames(df2) <- as.character(dest_priv$osm_id)
row.names(df2) <- as.character(ori$CODE_IRIS)
write.csv(df2, "data-conso/bike-duration-priv.csv")

# FFME
dest_asso <-dest[dest$type == "Associative structure",] 
dest_ffme <- dest_asso[dest_asso$federation == "FFME",]
df3 <- osrmTable(src = ori, dst = dest_ffme, measure = "duration")
df3 <- data.frame(df3$duration)
colnames(df3) <- as.character(dest_ffme$osm_id)
row.names(df3) <- as.character(ori$CODE_IRIS)
write.csv(df3, "data-conso/bike-duration-ffme.csv")

# FSGT
dest_fsgt <- dest_asso[dest_asso$federation == "FSGT",]
df4 <- osrmTable(src = ori, dst = dest_fsgt, measure = "duration")
df4 <- data.frame(df4$duration)
colnames(df4) <- as.character(dest_fsgt$osm_id)
row.names(df4) <- as.character(ori$CODE_IRIS)
write.csv(df4, "data-conso/bike-duration-fsgt.csv")

This heavy matrix is then manipulated to extract for each IRIS the following information :

  • Name of the nearest climbing structure by bike.

  • Type of manager (private, FFME or FSGT).

  • Time to reach it (in minutes).

  • Number of artificial climbing areas at less than 15 minutes by bike (N).

This calculation is done for all climbing structures (ALL), for private ones (PRIV), for those held by FFME federation (FFME) and for those held by the FSGT federation (FSGT). It is done with R base code. A function would reduce considerably the code length… Nevermind, it is not the aim of the proof ;)

Code 
# 5. Accessibility indicator creation (IRIS) ----
# Name of the nearest structure
df <- read.csv("data-conso/bike-duration.csv", row.names = "X")
colnames(df) <- as.character(dest$osm_id)
osm_id <- colnames(df)[apply(df, 1, which.min)] # Name
osm_id <- data.frame(osm_id, stringsAsFactors = FALSE)
osm_id$iris <- row.names(df)
osm_id <- merge(osm_id, poi[,c("osm_id", "name", "federation")], 
                by = "osm_id", all.x = TRUE)

# Time to the nearest climbing area
time <- apply(df, 1, min) # Time
time <- data.frame(time, stringsAsFactors = FALSE)
time$iris <- row.names(time)
osm_id <- merge(osm_id, time, by = "iris", all.x = TRUE)
osm_id$geometry <- NULL

# Number of climbing area at less than 15 minutes by bike
n15mn <- df
n15mn <- data.frame(df, stringsAsFactors = FALSE)
n15mn[n15mn <= 15] <- 1
n15mn[n15mn > 15] <- 0
n15mn$N <- rowSums(n15mn)
n15mn$iris <- row.names(n15mn)
osm_id <- merge(osm_id, n15mn[,c("iris", "N")], by = "iris",
                all.x = TRUE)
osm_id <- osm_id[,c(1,3:6)]
colnames(osm_id) <- c("CODE_IRIS", "ALL_NAME", "TYPE_STRUCT", "ALL_TIME",
                      "N_15MN")
iris <- merge(iris, osm_id, by = "CODE_IRIS", all.x = TRUE)

# Private climbing club (fees $$$)
# Name of the nearest structure
df2 <- read.csv("data-conso/bike-duration-priv.csv", row.names = "X")
colnames(df2) <- as.character(dest_priv$osm_id)
osm_id <- colnames(df2)[apply(df2, 1, which.min)] # Name
osm_id <- data.frame(osm_id, stringsAsFactors = FALSE)
osm_id$iris <- row.names(df2)
osm_id <- merge(osm_id, poi[,c("osm_id", "name", "type")], 
                by = "osm_id", all.x = TRUE)

# Time to the nearest climbing area
time <- apply(df2, 1, min) # Time
time <- data.frame(time, stringsAsFactors = FALSE)
time$iris <- row.names(time)
osm_id <- merge(osm_id, time, by = "iris", all.x = TRUE)
osm_id$geometry <- NULL

# Number of climbing area at less than 15 minutes by bike
n15mn <- df2
n15mn <- data.frame(df2, stringsAsFactors = FALSE)
n15mn[n15mn <= 15] <- 1
n15mn[n15mn > 15] <- 0
n15mn$N <- rowSums(n15mn)
n15mn$iris <- row.names(n15mn)
osm_id <- merge(osm_id, n15mn[,c("iris", "N")], by = "iris",
                all.x = TRUE)
osm_id <- osm_id[,c(1,3,5:6)]
colnames(osm_id) <- c("CODE_IRIS", "PRIV_NAME", "PRIV_TIME",
                      "N_PRIV_15MN")
iris <- merge(iris, osm_id, by = "CODE_IRIS", all.x = TRUE)

# FFME associative structure
# Name of the nearest structure
df3 <- read.csv("data-conso/bike-duration-ffme.csv", row.names = "X")
colnames(df3) <- as.character(dest_ffme$osm_id)
osm_id <- colnames(df3)[apply(df3, 1, which.min)] # Name
osm_id <- data.frame(osm_id, stringsAsFactors = FALSE)
osm_id$iris <- row.names(df3)
osm_id <- merge(osm_id, poi[,c("osm_id", "name", "type")], 
                by = "osm_id", all.x = TRUE)

# Time to the nearest climbing area
time <- apply(df3, 1, min) # Time
time <- data.frame(time, stringsAsFactors = FALSE)
time$iris <- row.names(time)
osm_id <- merge(osm_id, time, by = "iris", all.x = TRUE)
osm_id$geometry <- NULL

# Number of climbing area at less than 15 minutes by bike
n15mn <- df3
n15mn <- data.frame(df3, stringsAsFactors = FALSE)
n15mn[n15mn <= 15] <- 1
n15mn[n15mn > 15] <- 0
n15mn$N <- rowSums(n15mn)
n15mn$iris <- row.names(n15mn)
osm_id <- merge(osm_id, n15mn[,c("iris", "N")], by = "iris",
                all.x = TRUE)
osm_id <- osm_id[,c(1,3,5:6)]
colnames(osm_id) <- c("CODE_IRIS", "FFME_NAME", "FFME_TIME",
                      "N_FFME_15MN")
iris <- merge(iris, osm_id, by = "CODE_IRIS", all.x = TRUE)

# FSGT associative structure
# Name of the nearest structure
df4 <- read.csv("data-conso/bike-duration-fsgt.csv", row.names = "X")
colnames(df4) <- as.character(dest_fsgt$osm_id)
osm_id <- colnames(df4)[apply(df4, 1, which.min)] # Name
osm_id <- data.frame(osm_id, stringsAsFactors = FALSE)
osm_id$iris <- row.names(df4)
osm_id <- merge(osm_id, poi[,c("osm_id", "name", "type")], 
                by = "osm_id", all.x = TRUE)

# Time to the nearest climbing area
time <- apply(df4, 1, min) # Time
time <- data.frame(time, stringsAsFactors = FALSE)
time$iris <- row.names(time)
osm_id <- merge(osm_id, time, by = "iris", all.x = TRUE)
osm_id$geometry <- NULL

# Number of climbing area at less than 15 minutes by bike
n15mn <- df4
n15mn <- data.frame(df4, stringsAsFactors = FALSE)
n15mn[n15mn <= 15] <- 1
n15mn[n15mn > 15] <- 0
n15mn$N <- rowSums(n15mn)
n15mn$iris <- row.names(n15mn)
osm_id <- merge(osm_id, n15mn[,c("iris", "N")], by = "iris",
                all.x = TRUE)
osm_id <- osm_id[,c(1,3,5:6)]
colnames(osm_id) <- c("CODE_IRIS", "FSGT_NAME", "FSGT_TIME",
                      "N_FSGT_15MN")
iris <- merge(iris, osm_id, by = "CODE_IRIS", all.x = TRUE)

Characterizing the POI neighbourhood

Then, the socio-economic neighbourhood of each climbing area is characterized.

The previous matrix is transposed (lines <> columns). For each climbing structure, the IRIS located at less than 15 minutes by bike is extracted to produce the following indicators :

  • Total population at less than 15 minutes by bike. This can be considered as an amount of “local social demand”.

  • Minimum, mean and maximum of median income of IRIS at less than 15 minutes by bike.

Code 
# 6. Characterise the POI neighbourhood ----
t.df <- data.frame(t(df))
colnames(t.df) <- iris10k[iris10k$TYP_IRIS == "H",]$CODE_IRIS
t.df$osm_id <- dest$osm_id
t.df <- t.df[t.df$osm_id %in% poi$osm_id,]

poi_socio <- data.frame(matrix(nrow = 0, ncol = 5))
colnames(poi_socio) <- c("osm_id", "SUM_POP18", "MIN_REV19", "MOY_REV19", "MAX_REV19")

for (i in 1:nrow(t.df)){
  tmp1 <- t.df[, t.df[i, ] < 15]
  tmp2 <- data.frame(CODE_IRIS = colnames(tmp1))
  tmp2 <- merge(tmp2, iris10k[,c("CODE_IRIS", "P18_POP", "DISP_MED19")],
             all.x = TRUE)
  
  poi_socio[i,1] <- row.names(tmp1)[i]
  poi_socio[i,2] <- sum(tmp2$P18_POP, na.rm = TRUE)
  poi_socio[i,3] <- min(tmp2$DISP_MED19, na.rm = TRUE)
  poi_socio[i,4] <- mean(tmp2$DISP_MED19, na.rm = TRUE)
  poi_socio[i,5] <- max(tmp2$DISP_MED19, na.rm = TRUE)
  }

poi <- merge(poi, poi_socio, by = "osm_id", all.x = TRUE)

Output

Below is mapped one indicator calculated previously at IRIS level : time required to reach the nearest climbing area from the each populated IRIS centroid.

Code 
library(mapsf)

# Create the map
mf_init(com, expandBB = c(0,0.4,0,0), theme = my_theme)

cols <- mf_get_pal(n = 10, pal = "RdYlGn", rev = TRUE)
thr <- c(0, 2.5, 5, 7.5, 10, 12.5, 15, 20, 25, 30, max(iris$ALL_TIME, na.rm = TRUE))

mf_choro(
  x = iris, var = "ALL_TIME", 
  pal = cols, lwd = 0.25, breaks = thr, border = "white",
  leg_pos = c(631000, 6872000),
  leg_title = "Nearest climbing area\n(all structures, minutes)", 
  col_na = "lightgrey", leg_no_data = "No dedicated to housing", add = TRUE)

mf_map(x = poi, pch = 21, cex = .9, bg = "white", add = TRUE)
mf_map(com, col = NA, border = "black", add = TRUE)
mf_label(x = poi, var = "name", cex = .4, halo = TRUE, overlap = FALSE, 
         bg = "#ffffff80", pos = 4, add = TRUE)

mf_title("Time to reach the nearest climbing area")
mf_scale(size = 5, col = "black", pos = "bottomright")
mf_credits(paste0("Sources : © OpenStreetMap and Contributors, IGN, INSEE, 2022\n",
                  "Realisation : Ronan Ysebaert, 2022"))

Travel time isochrones

Moreover, we can be interested by defining the accessibility of climbing areas without taking into account any territorial division, which can introduce numerous bias, due to MAUP effects. This is done by building isochrones of equivalent time-distance coming from climbing areas, using the mapiso R package (Giraud (2022a)). Below are reminded the required steps to create polygons of equipotential from a regular grid of points :

  • Create a regular grid (150m resolution in our case) and extract the centroids.
  • Compute travel time indicators with the osrmtable function of the osrm R package, from grid centroids to climbing areas.
  • Get the minimum value of the resulting matrix and join the result to the regular grid layer.
  • Create polygons with mapiso function, wich takes the resulting values included in the regular grid. The function takes the desired thresholds and the indicator we are interested in (which climbing area travel time).
Code 
## 4.3 Isochrones from a regular grid to climbing areas ----
# Compute travel time from grid to climbing areas
# Create grid and extract centroids (cell size = 150 m)
mygrid <- st_make_grid(paris5k, cellsize = 150)
mygrid <- st_centroid(mygrid)
mygrid <- st_sf(ID = 1:length(mygrid), geometry = mygrid)

dest <- st_transform(dest, 2154)
dest_priv <- dest[dest$type == "Speculative structure",] 
dest_fsgt <- dest[dest$federation == "FSGT",]

# Compute travel time from grid centroids to all climbing areas
df5 <- osrmTable(src = mygrid, dst = dest, measure = "duration")
df5 <- data.frame(df5$duration)
write.csv(df5, "data-conso/grid-bike-duration.csv", row.names = FALSE)
df5 <- read.csv("data-conso/grid-bike-duration.csv")
colnames(df5) <- as.character(dest$osm_id)
time <- data.frame(mygrid$ID, apply(df5, 1, min)) # find minimum value
colnames(time) <- c("ID", "TIME_ALL")
mygrid <- merge(mygrid, time, by = "ID", all.x = TRUE) # merge time to grid

# Compute travel time from grid centroids to private climbing areas
df6 <- osrmTable(src = mygrid, dst = dest_priv, measure = "duration")
df6 <- data.frame(df6$duration)
write.csv(df6, "data-conso/grid-bike-duration_priv.csv", row.names = FALSE)
df6 <- read.csv("data-conso/grid-bike-duration_priv.csv")
colnames(df6) <- as.character(dest_priv$osm_id)
time <- data.frame(mygrid$ID, apply(df6, 1, min)) # find minimum value
colnames(time) <- c("ID", "TIME_PRIV")
mygrid <- merge(mygrid, time, by = "ID", all.x = TRUE) # merge time to grid

# Compute travel time from grid centroids to FSGT climbing areas
df7 <- osrmTable(src = mygrid, dst = dest_fsgt, measure = "duration")
df7 <- data.frame(df7$duration)
write.csv(df7, "data-conso/grid-bike-duration_fsgt.csv", row.names = FALSE)
df7 <- read.csv("data-conso/grid-bike-duration_fsgt.csv")
colnames(df7) <- as.character(dest_fsgt$osm_id)
time <- data.frame(mygrid$ID, apply(df7, 1, min)) # find minimum value
colnames(time) <- c("ID", "TIME_FSGT")
mygrid <- merge(mygrid, time, by = "ID", all.x = TRUE) # merge time to grid

# Compute isochrones
# define breaks (based on quantile analysis)
library(mapiso)
thr <- c(0, 2.5, 5, 7.5, 10, 12.5, 15, 20, 25, 30, max(mygrid$TIME_FSGT))
iso_all <- mapiso(x = mygrid, var = "TIME_ALL", breaks = thr, mask = paris5k)
iso_fsgt <- mapiso(x = mygrid, var = "TIME_FSGT", breaks = thr,  mask = paris5k)
iso_priv <- mapiso(x = mygrid, var = "TIME_PRIV", breaks = thr,  mask = paris5k)

Output

The maps below display the location of the climbing areas in the study area (white dots), the travel-time by bike required to reach each points from the 150m regular grid (dot map) and the resulting isochrones deducted from these values.

Code 
# thresholds, colours and background colours
thr <- c(0, 2.5, 5, 7.5, 10, 12.5, 15, 20, 25, 30, max(mygrid$TIME_ALL))

# Map input grid
mf_init(com, expandBB = c(0,0.4,0,0), theme = my_theme)
mf_map(mygrid, var = "TIME_ALL", pch = 21, cex = .4, border = NA,
       type = "choro", breaks = thr, pal = cols, leg_pos = c(631000, 6872000), 
       leg_title = "Time to reach\nthe nearest climbing area\nfrom a 150m regular\ngrid centroids",
       add = TRUE) 
mf_map(poi, pch = 21, add = TRUE)
mf_title("Isochrones input : travel time from a 150 m regular grid")
mf_scale(size = 5, col = "black", pos = "bottomright")
mf_credits(paste0("Sources : © OpenStreetMap and Contributors, 2022\n",
                  "Realisation : Ronan Ysebaert, 2022"))

Code 
# Map resulting isochrones
mf_init(com, expandBB = c(0,0.4,0,0), theme = my_theme)
mf_map(iso_all, var = "isomin", type = "choro", breaks = thr, pal = cols, 
       border = NA, leg_pos = c(631000, 6872000), 
       leg_title = "Time to reach\nthe nearest climbing area\nby bike",
       add = TRUE)
mf_map(com, col = NA, add = TRUE)
mf_map(poi, pch = 21, add = TRUE)
mf_title("Isochrones output : transformation of regularly spaced grids into contour polygons")
mf_scale(size = 5, col = "black", pos = "bottomright")
mf_credits(paste0("Sources : © OpenStreetMap and Contributors, IGN, 2022\n",
                  "Realisation : Ronan Ysebaert, 2022"))

OSRM climbing trip

osrmTrip is a function from the osrm package (Giraud (2022b)) which allows to get the travel geometry between multiple unordered points. The output is a sf LINESTRING with 2 components : duration (in minutes) and distance (in kilometres).

This function is applied to the climbing areas held by the FSGT federation.

Code 
dest_fsgt <- st_transform(dest_fsgt, 2154)
dest_fsgt <- st_intersection(dest_fsgt, paris_5k)
dest_fsgt <- st_transform(dest_fsgt, 4326)

trip <- osrmTrip(loc = dest_fsgt)
trip <- trip[[1]]$trip

To reach the 22 climbing areas of the FSGT federation located in the study area, the osrmTrip function returns that the fastest cycling trip around these climbing areas can be done in 420.93 minutes and 88.3846 kilometres.

Code 
library(maptiles)
Warning: le package 'maptiles' a été compilé avec la version R 4.2.1
Code 
osm <- get_tiles(x = trip, crop = TRUE, zoom = 13)

theme <- mf_theme(mar = c(0,0,1.2,0), inner = FALSE, line = 1.2, cex = .9, 
                  pos = "center", tab = FALSE)
mf_raster(osm)
mf_map(trip, lwd = 4, add = TRUE, col = "blue")
mf_map(trip, lwd = 1, col = "white", add = TRUE)
mf_map(dest_fsgt, pch = 20, col = "red", add = TRUE)
mf_title("Fastest cycling trip to reach all FSGT climbing areas")
mf_credits(get_credit("OpenStreetMap"), pos = "bottomright", 
           bg = "#ffffff80")

Manual corrections (out of OSM)

I noticed that in the OpenStreetMap database some attributes could be specified more precisely. However, my use does not correspond to the one recommended by OSM : climbing_length is not climbing_max. Moreover, this attribute has been set directly in OSM by the private brand.

Consequently, I prefer not to change these attributes in the OSM database, but directly in my R programme.

Code 
# Correct MurMur and Rename ESC15
poi[17,"climbing_length"] <- 17
poi[16,"name"] <- "ESC 15 - La Plaine"
poi[31,"name"] <- "ESC 15 - Croix Nivert"

Simplify geometries

Geometries are quite detailed. The library rmapshaper (Andy Teucher (2022)) is used to simplify the layer. 9 % of the points constituting the polygons are kept. Then, the IRIS layer is aggregated in communes for the map layout.

We also extract the IRIS located at less than 15 minutes from a artificial climbing area for some upcoming maps.

Code 
library(rmapshaper)
iris <- ms_simplify(iris, keep = 0.09)

# Communes aggregation (layout)
com <- aggregate(iris[,c("INSEE_COM", "NOM_COM")],
                 by = list(iris$INSEE_COM),
                 FUN = head, 1)

# Extract IRIS at less than 15 minutes by bike
iris15 <- iris[iris$ALL_TIME < 15,]
iris15 <- iris15[!is.na(iris15$ALL_TIME),]

Export results

Data preparation is over. Resulting sf objects that will be used for data visualizations are exported in geojson files in the data-conso folder.

Code 
# Export IRIS, POI, com layers
st_write(com, "data-conso/com.geojson")
st_write(iris, "data-conso/iris.geojson")
st_write(poi, "data-conso/poi.geojson")
st_write(iris15, "data-conso/iris15.geojson")

# Export material related to isochrones
mygrid <- st_transform(mygrid, 4326)
iso_all <- st_transform(iso_all, 4326)
iso_fsgt <- st_transform(iso_fsgt, 4326)
iso_priv <- st_transform(iso_priv, 4326)
st_write(mygrid, "data-conso/mygrid.geojson")
st_write(iso_all, "data-conso/iso_all.geojson")
st_write(iso_fsgt, "data-conso/iso_fsgt.geojson")
st_write(iso_priv, "data-conso/iso_priv.geojson")

# Export material related to travel-trip
st_write(trip, "data-conso/trip.geojson")

Session info (R)

Warning: le package 'osmdata' a été compilé avec la version R 4.2.1
Warning: le package 'readxl' a été compilé avec la version R 4.2.1
R version 4.2.0 (2022-04-22 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=French_France.utf8  LC_CTYPE=French_France.utf8   
[3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C                  
[5] LC_TIME=French_France.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] readxl_1.4.1     rmapshaper_0.4.6 osrm_4.0.0.0     osmdata_0.1.10  
[5] maptiles_0.3.0   mapsf_0.5.0      sf_1.0-8        

loaded via a namespace (and not attached):
 [1] tidyselect_1.1.2   terra_1.6-7        xfun_0.32          purrr_0.3.4       
 [5] lattice_0.20-45    V8_4.2.1           vctrs_0.4.1        generics_0.1.3    
 [9] htmltools_0.5.3    yaml_2.3.5         utf8_1.2.2         rlang_1.0.4       
[13] geojsonlint_0.4.0  e1071_1.7-11       pillar_1.8.1       glue_1.6.2        
[17] httpcode_0.3.0     DBI_1.1.3          slippymath_0.3.1   sp_1.5-0          
[21] lifecycle_1.0.1    stringr_1.4.1      cellranger_1.1.0   htmlwidgets_1.5.4 
[25] codetools_0.2-18   evaluate_0.16      knitr_1.40         fastmap_1.1.0     
[29] curl_4.3.2         class_7.3-20       fansi_1.0.3        Rcpp_1.0.9        
[33] KernSmooth_2.23-20 classInt_0.4-7     jsonvalidate_1.3.2 jsonlite_1.8.0    
[37] png_0.1-7          digest_0.6.29      stringi_1.7.8      dplyr_1.0.9       
[41] grid_4.2.0         cli_3.3.0          tools_4.2.0        magrittr_2.0.3    
[45] proxy_0.4-27       tibble_3.1.8       crul_1.2.0         pkgconfig_2.0.3   
[49] assertthat_0.2.1   rmarkdown_2.16     rstudioapi_0.14    R6_2.5.1          
[53] units_0.8-0        compiler_4.2.0

Data visualization with Observable JavaScript (ojs)

Observable is a startup founded by Mike Bostock and Melody Mechfessel, who initiated an on-line platform to collaboratively explore, analyse, visualize and communicate with data on the Web.

The computing language is Observable JavaScript (ojs). It is possible to include this programming language in Quarto chunks.

This section imports some visualizations realized in an Observable collection. Have a look to this resource to see the overall project !

All the maps produced below use the bertin.js library (Lambert (2023a)). Have a look to the documentation in this npm repository or in this Observable Collection to see the possibilities offered by this library to the map creator !

Import consolidated data

Geojson files coming from the data preparation precessing realized in a R programming language and described above are imported.

Code 
iris = FileAttachment("data-conso/iris.geojson").json()
poi = FileAttachment("data-conso/poi.geojson").json()
com = FileAttachment("data-conso/com.geojson").json()
iris15 = FileAttachment("data-conso/iris15.geojson").json()
iso_all = FileAttachment("data-conso/iso_all.geojson").json()
iso_priv = FileAttachment("data-conso/iso_priv.geojson").json()
iso_fsgt = FileAttachment("data-conso/iso_fsgt.geojson").json()
trip = FileAttachment("data-conso/trip.geojson").json()

Bertin.js (Lambert (2023a)), geotoolbox.js (Lambert (2023b)) geoverview.js (Lambert (2022)) libraries are imported. These libraries will be used for creating the maps.

Code 
bertin = require('bertin')
view = require("geoverview").then((f) => f.view)
geotoolbox = require('geotoolbox')

Layers attributes and metadata

Three geographical layers will be used :

  • iris : IRIS intersecting the study area (polygons layer).
  • poi : climbing structures included the study area (points layer).
  • com : Municipalities intersecting the study area (polygons layer). This layer will only be used for the map layout and better understand the location of the poi and iris layers : IRIS, the lowest territorial division in France, is not really known by the population.

The first two geographical layers include attributes coming from the data preparation process, that will be mapped and plotted in the next section of the document. In this part, we present the layers in several ways to understand the information we have in hand for producing the data visualization.

Interactive map

Made with geoverview library. Click on the points / polygons to have an overview of their respective attributes and values.

Code 
data = new Map([
  ["Climbing areas (dots)", poi],
  ["IRIS (polygons)", iris],
  ["Layout municipalities (polygons)", com]
])
viewof geojson = Inputs.select(data, { label: "Select a geojson" })
viewof style = Inputs.select(["voyager", "night", "fulldark","positron","icgc","osmbright","hibrid"], { label: "Select a style" })
map = view(geojson, {
  style: style
})

Attributes and metadata

Attribute table (iris layer)

Code 
Inputs.table(iris.features.map((d) => d.properties))

Indicators definition (iris layer)

id Description Data source
CODE_IRIS IRIS code INSEE, IGN
NOM_IRIS IRIS Name INSEE, IGN
TYP_IRIS IRIS type (H for dwellings) INSEE, IGN
P18_POP Population in 2018 INSEE
DISP_MED19 Disposable median income in 2019 INSEE
ALL_NAME Name of the nearest artificial climbing area reachable by bike (all structures) (c)OpenStreetMap and contributors, Ronan Ysebaert
TYPE_STRUCT Manager of the nearest artificial climbing area reachable by bike (private, FSGT or FFME) (c)OpenStreetMap and contributors, Ronan Ysebaert
ALL_TIME Time by bike to reach the nearest artificial climbing area. (c)OpenStreetMap and contributors, Ronan Ysebaert
N_15MN Number of climbing areas accessible in 15 minutes by bike from the IRIS centroid (all structures) (c)OpenStreetMap and contributors, Ronan Ysebaert
PRIV_NAME Name of the nearest artificial climbing area by bike (Private structure) (c)OpenStreetMap and contributors, Ronan Ysebaert
PRIV_TIME Time by bike to reach the nearest artificial climbing area (Private structure) (c)OpenStreetMap and contributors, Ronan Ysebaert
N_PRIV_15MN Number of climbing areas accessible in 15 minutes by bike from the IRIS centroid (private structures) (c)OpenStreetMap and contributors, Ronan Ysebaert
FFME_NAME Name of the nearest artificial climbing area by bike (Associative structure, federation : FFME) (c)OpenStreetMap and contributors, Ronan Ysebaert
FFME_TIME Time by bike to reach the nearest artificial climbing area (Associative structure, federation : FFME) (c)OpenStreetMap and contributors, Ronan Ysebaert
N_FFME_15MN Number of climbing areas accessible in 15 minutes by bike from the IRIS centroid (Associative structure, federation : FFME) (c)OpenStreetMap and contributors, Ronan Ysebaert
FSGT_NAME Name of the nearest artificial climbing area by bike (Associative structure, federation : FSGT) (c)OpenStreetMap and contributors, Ronan Ysebaert
FSGT_TIME Time by bike to reach the nearest artificial climbing area (Associative structure, federation : FSGT) (c)OpenStreetMap and contributors, Ronan Ysebaert
N_FSGT_15MN Number of climbing areas accessible in 15 minutes by bike from the IRIS centroid (Associative structure, federation : FSGT) (c)OpenStreetMap and contributors, Ronan Ysebaert

Attribute table (poi layer)

Code 
Inputs.table(poi.features.map((d) => d.properties))

Indicators definition (poi layer)

id Description Data source
osm_id OpenStreetMap identifier (responding to the query sport=climbing in the bounding box of the study area) (c)OpenStreetMap and contributors
name Name of the artificial climbing area (c)OpenStreetMap and contributors
climbing_type If the artificial climbing area is a wall (need a rope to climb), dedicated to bouldering, or both (c)OpenStreetMap and contributors, Ronan Ysebaert
climbing_length Max. height of the wall (if adapted, for walls and not for bouldering) (c)OpenStreetMap and contributors
climbing_routes Number of climbing routes, materialized by belay station (if adapted, for walls and not for bouldering). Note that in artificial climbing, a same belay station can serve several climbing routes (c)OpenStreetMap and contributors
type If the climbing area is speculative (manage by a private brand) or associative (c)OpenStreetMap and contributors, Ronan Ysebaert
federation Name of the manager (FFME, FSGT or private) (c)OpenStreetMap and contributors, Ronan Ysebaert
brand Name of the brand managing the private structure (c)OpenStreetMap and contributors
SUM_POP18 Population at less than 15 minutes by bike of the artificial climbing area (IRIS centroids). INSEE, (c)OpenStreetMap and contributors, Ronan Ysebaert
MIN_REV19 Minimum of the median income for the IRIS located at less than 15 minutes by bike INSEE, (c)OpenStreetMap and contributors, Ronan Ysebaert
MOY_REV19 Mean of the Median income for the IRIS located at less than 15 minutes by bike INSEE, (c)OpenStreetMap and contributors, Ronan Ysebaert
MAX_REV19 Maximum of the Median income for the IRIS located at less than 15 minutes by bike INSEE, (c)OpenStreetMap and contributors, Ronan Ysebaert

Wall’s characteristics

This interactive map displays the artificial climbing offer in Paris and its surroundings, according artificial climbing specificities: top rope climbing (walls) / bouldering, climbing height, number of routes, manager of the artificial climbing area.

Code 
lambert93 = "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

// Map parameters to play with it
viewof ind = Inputs.radio(["climbing_length", "climbing_routes"], {label: "Wall height / Number of routes", value: "climbing_length"})
viewof k = Inputs.range([10, 300], { label: "Spikes height", step: 1, value: 90 })
viewof w = Inputs.range([2, 30], { label: "Spikes / dots width", step: 1, value: 15 })

// Spike map with tooltips
map1 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },

  layers: [
    {
      type: "spikes",
      geojson: poi,
      values: ind,
      fillOpacity: 0.2,
      k: k,
      w: w,
      fill: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"],
        leg_title: `Wall's managers`,
        leg_x: 600,
        leg_y: 550
      },
      stroke: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"]
      },
      leg_x: 15,
      leg_y: 15,
      leg_round: 0,
      leg_title: `Height (meters) or routes number`,
      tooltip: {
        fields: [
          "$name",
          "$climbing_length",
          "(height in meters)",
          "$climbing_routes",
          "(nb. routes)",
          "$federation",
          "(manager)"
        ],
        fill: ["black", "black", "black", "black", "black", "black", "black"],
        fontWeight: [
          "bold",
          "normal",
          "normal",
          "normal",
          "normal",
          "normal",
          "normal"
        ],
        fontSize: [14, 12, 10, 12, 10, 12, 10],
        col: "white"
      }
    },

    {
      geojson: poi,
      symbol_size: 10 * w,
      symbol: "circle",
      fill: {
        type: "typo",
        values: "climbing_type",
        colors: ["white", "lightgrey", "grey"],
        leg_title: `Climbing type`,
        leg_x: 600,
        leg_y: 430
      },
      stroke: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"]
      },
      strokeWidth: 2,
      tooltip: {
        fields: ["$name", "$climbing_type"],
        fill: ["black", "black"],
        fontWeight: ["bold", "normal"],
        fontSize: [12, 10],
        col: "white"
      }
    },

    {
      type: "text",
      position: "topright",
      text: "Climbing hot spots",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "black",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      type: "layer",
      geojson: com,
      stroke: "white",
      fill: "peachpuff",
      tooltip: ["$NOM_COM"]
    },

    {
      type: "scalebar",
      units: "miles"
    },

    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Nearest wall

The nearest climbing area by bike from each populated IRIS by wall manager: private or associative (FSGT or FFME).

Code 
// Categorical map with tooltips
map2 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },

  layers: [
    {
      geojson: poi,
      symbol_size: 100,
      symbol: "circle",
      fill: "white",
      strokeWidth: 2.5,
      stroke: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"]
      },
      
      tooltip: {
        fields: ["$name", "$federation"],
        fill: ["black", "black"],
        fontWeight: ["bold", "normal"],
        fontSize: [12, 10],
        col: "white"
      }
    },
      
    {
      type: "text",
      position: "topright",
      text: "Nearest climbing area\nmanager",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "black",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      type: "layer",
      geojson: com,
      stroke: "black",
      fill: "none"
    },

    {
      type: "layer",
      geojson: iris,
      fill: {
        type: "typo",
        values: "TYPE_STRUCT",
        strokeWidth: 0.3,
        colors: ["#fcb45680", "#ba0b2580", "#9990f080"],
        leg_title: `Manager name`,
        leg_x: 15,
        leg_y: 15,
        col_missing: "lightgrey",
        txt_missing: "Non dedicated to housing"
      },
      tooltip: {
        fields: [
          "$NOM_IRIS",
          "$ALL_NAME",
          "Name of the nearest artificial climbing structure"
        ],
        fill: "white",
        fontWeight: ["bold", "normal", "normal"],
        fontSize: [11, 10, 8],
        col: ["black", "black", "black"]
      }
    },

    {
      type: "scalebar",
      units: "miles",
    },
    
    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Time to climb (IRIS)

It displays the time by bike required to reach the nearest climbing area, for all structures or by manager type and from each IRIS centroids.

Code 
viewof ind2 = Inputs.radio(["ALL_TIME", "PRIV_TIME", "FSGT_TIME", "FFME_TIME"],
{value: "ALL_TIME"})

// Discretization method
viewof method2 = Inputs.select(["jenks", "q6", "quantile", "equal", "msd"], {
  label: "method",
  value: "quantile"
})

// nbreaks
viewof nbreaks2 = Inputs.range([3, 9], { label: "nbreaks", step: 1, value: 8 })

// colour palettes
viewof pal2 = Inputs.select(
  ["Blues",  "Greens",  "Greys", "Oranges", "Purples", "BuGn", "BuPu", "GnBu",
  "OrRd", "PuBuGn", "PuBu", "PuRd", "RdPu", "YlGnBu", "YlGn", "YlOrBr", "YlOrRd",
  "BrBG", "PRGn", "PiYG", "PuOr", "RdBu", "RdGy", "RdYlBu", "RdYlGn", "Spectral" ],
  {
    label: "colors",
    value: "RdYlBu"
  }
)

// Choropleth map with tooltips
map3 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },

  layers: [
    {
      geojson: poi,
      symbol_size: 100,
      symbol: "circle",
      fill: "white",
      strokeWidth: 2.5,
      stroke: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"]
      },
      tooltip: {
        fields: ["$name", "$federation"],
        fill: ["black", "black"],
        fontWeight: ["bold", "normal"],
        fontSize: [12, 10],
        col: "white"
      }
    },

    {
      type: "text",
      position: "topright",
      text: "Time to reach the nearest\nartificial climbing area",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "black",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      type: "layer",
      geojson: com,
      stroke: "black",
      fill: "none"
    },

    {
      type: "layer",
      geojson: iris,
      fill: {
        type: "choro",
        values: ind2,
        nbreaks: 8,
        method: method2,
        colors: pal2,
        leg_round: 2,
        leg_title: `By bike (minutes)`,
        leg_x: 15,
        leg_y: 15,
        col_missing: "lightgrey",
        txt_missing: "Non dedicated to housing"
      },
      tooltip: {
        fields: [
          "$NOM_IRIS",
          "$NOM_COM",
          "$ALL_TIME",
          "Bike time to the nearest Indoor climbing area"
        ],
        fill: "white",
        fontWeight: ["bold", "normal", "normal", "normal"],
        fontSize: [11, 11, 10, 8],
        col: ["black", "black", "black", "black"]
      }
    },

    {
      type: "scalebar",
      units: "kilometers",
    },
    
    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Time to reach climbing areas (isochrones)

It shows with isochrones the area reachable from the network of climbing areas by bike in minutes. It is possible to see all climbing areas (default map) or filter by manager (private or FSGT).

Code 
viewof ind7 = Inputs.radio([iso_all, iso_priv, iso_fsgt], {format: x => x.name, value: iso_all})

map7 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },

  layers: [
    {
      type: "layer",
      geojson: poi,
      symbol_size: 100,
      symbol: "circle",
      fill: "white",
      strokeWidth: 2.5,
      stroke: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"]
      },
      tooltip: {
        fields: ["$name", "$federation"],
        fill: ["black", "black"],
        fontWeight: ["bold", "normal"],
        fontSize: [12, 10],
        col: "white"
      }
    },
    {
      type: "text",
      position: "topright",
      text: "Climbing nebula",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "black",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      type: "layer",
      geojson: com,
      stroke: "black",
      fill: "none"
    },

    {
      type: "layer",
      geojson: ind7,
      strokeWidth: 0,
      fill: {
        type: "choro",
        values: "isomin",
        breaks: [0, 2.5, 5, 7.5, 10, 12.5, 15, 20, 25, 30, 47],
        colors: ["#006228", "#3A8B3A", "#76B345", "#B4D66C", "#EAF3A3",
        "#FAEDA9", "#F1C363", "#E98E00", "#DC4A00", "#A51122"],
        leg_round: 1,
        leg_title: `Time to reach\nthe nearest climbing area\nBy bike (minutes)`,
        leg_x: 15,
        leg_y: 15,
        col_missing: "lightgrey",
      }
    },

    {
      type: "scalebar",
      units: "kilometers"
    },

    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Climbing in 15 minutes?

The number of climbing areas reachable in 15 minutes by bike from the geometric centroid of each IRIS of the study area.

Code 
viewof ind4 = Inputs.radio(["N_15MN", "N_PRIV_15MN", "N_FSGT_15MN", "N_FFME_15MN"],
{label: "Number - all / private / FSGT / FFME - indoor climbing areas at less than 15 minutes", value: "N_15MN"})

// colour palettes
viewof pal4 = Inputs.select(
  ["Blues", "Greens", "Greys", "Oranges", "Purples", "BuGn", "BuPu", "GnBu", 
  "OrRd", "PuBuGn", "PuBu", "PuRd", "RdPu",
  ],
  {
    label: "colors",
    value: "Greens"
  }
)

// choropleth map with tooltips
map4 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },

  layers: [
    {
      geojson: poi,
      symbol_size: 100,
      symbol: "circle",
      fill: "white",
      strokeWidth: 2.5,
      stroke: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"]
      },
      
      tooltip: {
        fields: ["$name", "$federation"],
        fill: ["black", "black"],
        fontWeight: ["bold", "normal"],
        fontSize: [12, 10],
        col: "white"
      }
    },
      
    {
      type: "text",
      position: "topright",
      text: "Number of climbing areas reachable\nin less than 15 minutes",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "black",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      type: "layer",
      geojson: com,
      stroke: "black",
      fill: "none"
    },

    {
      type: "layer",
      geojson: iris,
      fill: {
        type: "choro",
        values: ind4,
        breaks: [1, 2, 3, 4, 5, 6, 7, 8, 9],
        colors: pal4,
        leg_round: 2,
        leg_title: `Structures numbers`,
        leg_x: 15,
        leg_y: 15,
        col_missing: "lightgrey",
        txt_missing: "Non dedicated to housing"
      },
      tooltip: {
        fields: [
          "$NOM_IRIS",
          "$N_15MN",
          "Climbing areas at less than 15 minutes by bike"
        ],
        fill: "white",
        fontWeight: ["bold", "normal", "normal"],
        fontSize: [11, 10, 8],
        col: ["black", "black", "black"]
      }
    },

    {
      type: "scalebar",
      units: "miles",
    },
    
    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Climbing / biking tour !

It shows the fastest path to visit all the FSGT’s climbing areas by bike. To overcome this challenge, you will have to cycle for 88.39 km and 7.02 hours !

Code 
viewof ind8 = Inputs.radio(["duration", "distance"], {
  label: "Choose bike trip indicator",
  value: "duration"
})

// Discretization method
viewof method8 = Inputs.select(["jenks", "q6", "quantile", "equal", "msd"], {
  label: "method",
  value: "equal"
})

// Filter FSGT climbing areas only
fsgt = geotoolbox.filter(poi, (d) => d.federation == "FSGT")


//fsgt = bertin.subgeo({
//  geojson: poi,
//  field: "federation",
//  operator: "==",
//  array: "FSGT"
//})

// Colour palette
viewof pal8 = Inputs.select(
  ["Reds", "Blues", "Greens", "Greys", "Oranges", "Purples", "BuGn", "BuPu", "GnBu",
  "OrRd", "PuBuGn", "PuBu", "PuRd", "RdPu", "YlGnBu", "YlGn", "YlOrBr", 
  "YlOrRd", "BrBG", "PRGn", "PiYG", "PuOr", "RdBu", "RdGy", "RdYlBu", "RdYlGn",
  "Spectral"
  ],
  {
    label: "colors",
    value: "Reds"
  }
)


// Resulting map with a lot of things :) 
map8 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },
  layers: [
    {
      geojson: fsgt,
      symbol_size: 100,
      symbol: "circle",
      fill: "white",
      strokeWidth: 2.5,
      stroke: "#ba0b25",

      tooltip: {
        fields: ["$name", "$federation"],
        fill: ["black", "black"],
        fontWeight: ["bold", "normal"],
        fontSize: [12, 10],
        col: "white"
      }
    },

    {
      type: "text",
      position: "topright",
      text: "The FSGT's\nbiking / climbing tour !",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "white",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      geojson: trip,
      strokeWidth: 4,
      stroke: {
        type: "choro",
        values: ind8,
        method: method8,
        nbreaks: 4,
        colors: pal8,
        leg_title: `Time (minutes) or distance (kilometres)\nbetween 2 climbing areas !`,
        leg_x: 15,
        leg_y: 15,
        leg_txtcol: "white",
        leg_round:0
      },

      tooltip: {
        fields: ["$duration","minutes", "$distance", "kilometers"],
        fill: ["black", "black"],
        fontWeight: ["bold", "normal"],
        fontSize: [12, 10, 12, 10],
        col: "white"
      }
    },

    {
      type: "layer",
      geojson: com,
      stroke: "white",
      fill: "black"
    },
    
    {
      type: "scalebar",
      units: "kilometers",
      fill: "white"
    },

    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Demographic neighbourhood

Visualize the demographic neighbourhood of each climbing area : IRIS population at less than 15 minutes by bike from each climbing area.

Code 
viewof k5 = Inputs.range([1, 100], { label: "radius max", step: 1, value: 18 })

// Proportional map and tooltips
map5 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },

  layers: [
    {
      geojson: poi,
      type: "bubble",
      k: k5,
      values: "SUM_POP18",
      leg_title: `Population in 15mn by bike neighbourhood`, 
      leg_x: 15,
      leg_y: 15,
      leg_round: -2,
      fill: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"],
        leg_title: `Wall's managers`,
        leg_x: 600,
        leg_y: 550
      },
      strokeWidth: 2.5,
      stroke: "white",
      tooltip: {
        fields: ["$name", "$federation", "$SUM_POP18", "(Inhabitants in 15mn by bike neighbourhood)"],
        fill: ["black", "black", "black", "black"],
        fontWeight: ["bold", "normal", "normal", "normal"],
        fontSize: [14, 12, 12, 10],
        col: "white"
      }
    },

    {
      type: "text",
      position: "topright",
      text: "Population\nin climbing area's neighbourhood",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "black",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      type: "layer",
      geojson: com,
      stroke: "black",
      fill: "none"
    },

    {
      type: "layer",
      geojson: iris15,
      fill: "peachpuff",
      tooltip: {
        fields: [
          "$NOM_IRIS",
          "$NOM_COM",
          "$P18_POP",
          "Total population (2018)"
        ],
        fill: "white",
        fontWeight: ["bold", "normal", "normal", "normal"],
        fontSize: [11, 11, 10, 8],
        col: ["black", "black", "black", "black"]
      }
    },

    {
      type: "scalebar",
      units: "kilometers",
    },
    
    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Income neighbourhood

Visualize the level of wealth of each climbing area neighbourhood : IRIS median income at less than 15 minutes by bike from a climbing area. User can choose the minimum income, average income or maximum income of the IRIS at less than 15 minutes.

Code 
viewof ind6 = Inputs.radio(["MIN_REV19", "MOY_REV19", "MAX_REV19"], {label: "Min/average/max median income in a 15 minutes by bike neighbourhood", value: "MOY_REV19"})

// Discretization method
viewof method6 = Inputs.select(["jenks", "q6", "quantile", "equal", "msd"], {
  label: "method",
  value: "quantile"
})

// Nbreaks
viewof nbreaks6 = Inputs.range([3, 9], { label: "nbreaks", step: 1, value: 6 })

// Colour palette
viewof pal6 = Inputs.select(
  ["Blues", "Greens", "Greys", "Oranges", "Purples", "BuGn", "BuPu", "GnBu",
  "OrRd", "PuBuGn", "PuBu", "PuRd", "RdPu", "YlGnBu", "YlGn", "YlOrBr", 
  "YlOrRd", "BrBG", "PRGn", "PiYG", "PuOr", "RdBu", "RdGy", "RdYlBu", "RdYlGn",
  "Spectral"
  ],
  {
    label: "colors",
    value: "PuOr"
  }
)

// IRIS layer opacity
viewof opacity6 = Inputs.range([0, 1], {label: "IRIS income opacity", step: 0.1})


// Resulting map with a lot of things :) 
map6 = bertin.draw({
  params: {
    width: 750,
    projection: lambert93,
    extent: com
  },

  layers: [
    {
      geojson: poi,
      symbol_size: 300,
      symbol: "triangle",
      fill: {
        type: "choro",
        values: ind6,
        nbreaks: nbreaks6,
        method: method6,
        colors: pal6,
      },
      strokeWidth: 2.5,
      stroke: {
        type: "typo",
        values: "federation",
        colors: ["#9990f0", "#fcb456", "#ba0b25"],
        leg_title: `Wall's managers`,
        leg_x: 600,
        leg_y: 550
      },
      tooltip: {
        fields: ["$name", "$federation", "$MIN_REV19", "(Minimum income in 15mn by bike neighbourhood)", "$MOY_REV19", "(Average income in 15mn by bike neighbourhood)", "$MAX_REV19", "(Average income in 15mn by bike neighbourhood)"],
        fill: ["black", "black", "black", "black", "black", "black", "black", "black"],
        fontWeight: ["bold", "normal", "normal", "normal", "normal", "normal", "normal", "normal"],
        fontSize: [14, 12, 12, 10, 12, 10, 12, 10],
        col: "white"
      }
    },

    {
      type: "text",
      position: "topright",
      text: "Income\nin climbing area's neighbourhood",
      frame_opacity: 0.7,
      baseline: "hanging",
      fill: "black",
      fontSize: 25,
      margin: 4,
      fontWeight: "bold",
      fontFamily: "Ubuntu"
    },

    {
      type: "layer",
      geojson: com,
      stroke: "black",
      fill: "none"
    },

    {
      type: "layer",
      geojson: iris15,
      fill: {
        type: "choro",
        values: "DISP_MED19",
        nbreaks: nbreaks6,
        method: method6,
        colors: pal6,
        leg_round: -2,
        leg_title: `Median income (2019, euros)\nfor IRIS at less than 15km by bike\nfrom a climbing area`,
        leg_x: 15,
        leg_y: 15,
        col_missing: "lightgrey",
        txt_missing: "Non dedicated to housing"
      },
      fillOpacity: opacity6,
      tooltip: {
        fields: [
          "$NOM_IRIS",
          "$NOM_COM",
          "$DISP_MED19",
          "Median income (2019, euros)"
        ],
        fill: "white",
        fontWeight: ["bold", "normal", "normal", "normal"],
        fontSize: [11, 11, 10, 8],
        col: ["black", "black", "black", "black"]
      }
    },

    {
      type: "scalebar",
      units: "kilometers",
    },
    
    {
      type: "footer",
      text:
        "Sources: (c)OpenStreetMap and contributors, INSEE, IGN, 2022\nRealisation: Ronan Ysebaert, 2022"
    }
  ]
})

Concluding methodological remarks

Extend the analysis

This methodological example showed a general methodology to combine institutional and crowd-sourced data to develop some simple indicators of accessibility to services of interest using open-source routing engines, and disseminate this with data visualizations frameworks. It is one possible implementation, among others. Other implementations could easily be envisaged, such as:

  • Geographical coverage: The data preparation process presented above requires the definition of a geographical bounding box to extract OpenStreetMap points of interest (destinations of travel-time calculations) and an existing territorial division (the centroids being the origins of the travel-time calculations). This could be easily extended to other areas of interest, by modifying the Xmin, Xmax, Ymin and Ymax of the bounding box, or using other territorial divisions (a regular grid, for instance). The important things to consider is
  • Use other points of interest: This example used the OpenStreetMap tag defined by the key-value sport=climbing. Other map features could be investigated, such as amenities, tourism facilities or shops. The only element to consider is the geometric reference of the object in the OpenStreetMap database ; it could be nodes, ways or relations. In the case of sport=climbing tag, the OpenStreetMap Wiki reminds that sport=climbing should preferably applied to nodes for artificial climbing walls, our case. If the tag combines nodes and areas, it requires also to extract polygons responding to the tag and transform it into points, checking if the information is not included twice in the OpenStreetMap database.
  • Apply other profiles: In this example, the standard bike profile has been applied. The pedestrian and the car profile are also available and may be also used depending on the objective of the analysis. Whatever the profile used, it reminds important to the parameters of the profiles (car here, bike here and pedestrian here) to understand the calculations realized with routing engines, and their limitations. It is also possible to implement its own profile.

Usage precautions

One persisting criticism made to OpenStreetMap is its lack of quality, mainly due to its collaborative aspect implying de facto a high heterogeneity of contribution practices. Known problems refer to geometric accuracy (geographical positioning of OSM objects), accuracy, completeness (missing points, not up-to-date) or problems due to the tags describing OSM objects (missing important tags, bad practices).

In that case, I made a personal survey to complete missing information (around 30 % of the information was missing) in OpenStreetMap database. It meant add some missing points in some cases, and in others to describe it with relevant climbing tags identified in the OpenStreetMap Wiki. At the end, it is sure that the input information is not perfect and would require double-checks by the community, but we can argue that it was not worst than in the beginning.

Whatever the use of OpenStreetMap you may made, taking a step back to these questions is very important to control a minima the input information used for the travel-time calculations.

Another point to have in mind is the default parameters of osrm profiles: bike profile does not take into account the slope of the road ; the car profile does not take into account the traffic congestion. It may have significant impact on the travel-time calculations.

Nevertheless, we argue that OpenStreetMap is a wonderful initiative and as we try to demonstrate in this notebook, it is not too difficult to produce innovative indicators and data visualizations using open-source database and libraries.

Enjoy !

References

Andy Teucher, Kenton Russell. 2022. Rmapshaper: Client for ’Mapshaper’ for ’Geospatial’ Operations. https://cran.r-project.org/web/packages/rmapshaper/index.html.
Giraud, Timothée. 2022a. Mapiso: Create Contour Polygons from Regular Grids. https://cran.r-project.org/web/packages/mapiso/index.html.
———. 2022b. Osrm: Interface Between r and the OpenStreetMap-Based Routing Service OSRM. https://CRAN.R-project.org/package=osrm.
Lambert, Nicolas. 2022. Geoverview.js: A Tool for Giving a Quick and Easy Geographic Overview of Any Geojson. https://www.npmjs.com/package/geoverview.
———. 2023a. Bertin.js: A JavaScript Library for Visualizing Geospatial Data and Make Thematic Maps for the Web. https://www.npmjs.com/package/bertin.
———. 2023b. Geotoolbox.js: A Tool for Geographers to Manage Properties (Attribute Data) and Provide Several Useful GIS Operations for Thematic Cartography. https://www.npmjs.com/package/geotoolbox.
Mark Padgham, Robin Lovelace, Bob Rudis. 2022. Osmdata: Import ‘OpenStreetMap‘ Data as Simple Features of Spatial Objects. https://CRAN.R-project.org/package=osmdata.

Citation

BibTeX citation:
@online{ysebaert,
  author = {Ronan Ysebaert},
  editor = {},
  title = {Travel Time Calculation with {R} and Data Visualization with
    {Observable}},
  date = {},
  langid = {en}
}
For attribution, please cite this work as:
Ronan Ysebaert. n.d. “Travel Time Calculation with R and Data Visualization with Observable.”