Railways

library(sf)library(glue)library(ggplot2)library(janitor)library(ggspatial)library(dplyr)library(tidyr)library(units)library(osmdata)

# Get and cache OSM data for abandoned railways in Franceif (!file.exists("lost_railways.gpkg")) {  railways <- r"([out:xml][timeout:6000];    // Get France as area  area["name"="France"]->.fr;    (  // Direct generic railway status tags  way["railway"~"^(disused|abandoned|razed|dismantled|demolished|retired)$"](area.fr);    // Include specific rail / tram / light_rail... in same statuses  way["railway"~"^(rail|tram|light_rail|narrow_gauge)$"]["disused"="yes"](area.fr);  way["railway"~"^(rail|tram|light_rail|narrow_gauge)$"]["abandoned"="yes"](area.fr);  way["railway"~"^(rail|tram|light_rail|narrow_gauge)$"]["razed"="yes"](area.fr);  way["railway"~"^(rail|tram|light_rail|narrow_gauge)$"]["dismantled"="yes"](area.fr);  way["railway"~"^(rail|tram|light_rail|narrow_gauge)$"]["demolished"="yes"](area.fr);  way["railway"~"^(rail|tram|light_rail|narrow_gauge)$"]["retired"="yes"](area.fr);    // Lifecycle prefix tags for some railway subtype  // Matches patterns like: disused:railway=rail, abandoned:railway=tram, etc.  way[~"^(disused|abandoned|razed|dismantled|demolished|retired):railway$"~"^(rail|tram|light_rail|narrow_gauge)$"](area.fr);  );  (._;>;);  out body;)" |>     osmdata_sf()     lost_railways <- railways$osm_lines |>     clean_names() |>     select(railway, name) |>     write_sf("lost_railways.gpkg")} else {  lost_railways <- read_sf("lost_railways.gpkg")}# Hexagonal grid# See https://r.iresmi.net/posts/2024/discrete_global_grid/hex_size <- 20 # kmhex <- read_sf("../../2024/discrete_global_grid/dggrid_fx.gpkg",                 layer = glue("dggrid_{hex_size}k")) |>   st_transform("EPSG:2154")dep <- "~/data/adminexpress/adminexpress_cog_simpl_000_2022.gpkg" |>   read_sf(layer = "departement_int")

railways_length <- lost_railways |>   st_transform("EPSG:2154") |>   st_intersection(hex) |>   mutate(length = st_length(geom)) |>   st_drop_geometry() |>   group_by(seqnum) |>   summarise(length = drop_units(set_units(sum(length), km)))

hex |>   left_join(railways_length,            join_by(seqnum)) |>  replace_na(list(length = 0)) |>   ggplot() +  geom_sf(aes(fill = length, color = length)) +  geom_sf(data = dep, fill = NA, color = "lightgrey", linewidth = 0.3, alpha = 0.8) +  scale_fill_viridis_c(name = "former railways length\nper cell (km)",                        trans = "sqrt",                       aesthetics = c("fill", "color")) +  annotation_scale(height = unit(1, "mm"),                    text_col = "darkgrey", line_col = "grey",                    bar_cols = c("white", "grey")) +  labs(title = "The vanished railways",       subtitle = glue("Currently disused, abandoned, dismantled... \\                       railways in France"),       caption = glue("data: OpenStreetMap contributors                      départements based on IGN AdminExpress 2022                      Discrete Global Grid ISEA3H ≈ {hex_size} km,                       one hexagon ≈ {round(hex_size^2 * sqrt(3) / 2)} km²                      https://r.iresmi.net - {Sys.Date()}")) +  theme_void() +  theme(plot.caption = element_text(size = 7, color = "grey40"),        plot.margin = unit(c(.2, .2, .2, .2), units = "cm"),        legend.position = "bottom")

Figure 1: Length per hexagon cell of currently disused, abandoned, dismantled, etc. railways in France