Skip to content

Instantly share code, notes, and snippets.

@szechno
Created February 6, 2026 11:56
Show Gist options
  • Select an option

  • Save szechno/d22428d956bc99883042e740ef30f067 to your computer and use it in GitHub Desktop.

Select an option

Save szechno/d22428d956bc99883042e740ef30f067 to your computer and use it in GitHub Desktop.
# prelim functions
# download oemextract and clean up ----
get_oemextract <- function(wsx = "D:/data_store/Boundaries/county/county.shp"){
message("Warning this will take 5-10 mins")
wsx <- sf::read_sf(wsx) |> dplyr::select(geometry)
wsx_roads <- osmextract::oe_get(
wsx |> st_transform(crs = 4326),
query = "SELECT * FROM 'lines' WHERE highway IN ('footway', 'unclassified',
'service', 'track', 'bridleway', 'path', 'residential', 'tertiary',
'primary', 'trunk', 'secondary', 'motorway', 'tertiary_link', 'pedestrian',
'cycleway', 'primary_link', 'motorway_link', 'trunk_link', 'secondary_link',
'living_street', 'construction', 'services', 'proposed', 'busway',
'crossing')"
)
message(paste("Continuing with intersection...", Sys.time()))
a <- Sys.time()
wsx_roads <- wsx_roads |>
sf::st_intersection(wsx |> st_transform(crs = 4326))
b <- Sys.time()
message(paste("sf::st_intersection:", b - a))
# "sf::st_intersection: 5.3903" as 4326 and wsx 27700
# "sf::st_intersection: 4.48411624828974" as both 4326.. fastest
# and as both 27700...
# didn't complete after 10 mins
# neither did st_intersection_faster and geos_intersection
wsx_roads <- wsx_roads |> st_transform(crs = 27700)
if(any(st_geometry_type(wsx_roads) != "LINESTRING")){
message("Converting multilinestrings to linestrings...")
wsx_roads <- wsx_roads |>
st_cast("MULTILINESTRING") |>
st_cast("LINESTRING")
}
return(wsx_roads)
}
# make wsx_roads sfnetwork ----
make_sfnetwork <- function(roads = wsx_roads){
message("creating sfnetwork...")
roads_sfn = as_sfnetwork(wsx_roads, directed = FALSE)
message("cleaning up sfnetwork...")
roads_clean <- roads_sfn |>
convert(to_spatial_subdivision, .clean = TRUE) |>
convert(to_spatial_smooth, .clean = TRUE)
return(roads_clean)
}
# start helper functions ----
# 0. list all unique schools in pupils ----
schools_list <- function(targets = pupils){
targets$nativeid |> unique() |> sort()
}
# 1. filter source ----
select_source <- function(source_id = "1000",
source = schools){
source |>
filter(nativeid == source_id) |>
mutate(
type = "source"
)
}
# 2. buffer network around source ----
buffer_roads <- function(source,
roads = roads_clean,
dist = BUFFER_DIST){
if(nrow(source) > 1 | is.null(source)){
warning("this source not 1 or is null")
return(NULL)
} else {
roads |> st_intersection(
source |>
select(geometry) |>
st_buffer(
dist = dist
)
)
}
}
# 3. filter targets by source ----
select_all_targets <- function(source,
all_targets = pupils,
dist = WALKING_DIST){
targets <- all_targets |>
filter(nativeid == source$nativeid) |>
mutate(
# UID = row_number(),
type = "target",
Count = 1
) |> st_intersection(
source |>
select(geometry) |>
st_buffer(
dist = dist
)
)
}
# 4. blend this school to buffered network ----
blend_points_to_sfnetwork <- function(source,
targets,
roads = roads_buff){
points <- bind_rows(source, targets)
road <- st_network_blend(roads, points, tolerance = Inf) |>
activate(nodes)|>
mutate(id = 1:n())
}
# 5. pull ids to use ----
# from is the source node id
# to is each of the target node ids
# ccount is used to carry over the count for each edge in the route
pull_ids <- function(roads = school_roads){
from <- roads |> activate(nodes) |> filter(type == "source") |> pull(id)
to <- roads |> activate(nodes) |> filter(type == "target") |> pull(id)
ccount <- roads |> activate(nodes) |> filter(type == "target") |> pull(Count)
return(list(from = from, to = to, ccount = ccount))
}
# 6. get paths ----
# used by get_routes
get_paths <- function(i,
roads = school_roads,
ids = pulled_ids){
paths <- st_network_paths(roads, from = ids$from, to = ids$to[[i]])
route <- roads |>
activate("edges")|>
slice(unlist(paths$edge_paths)) |>
st_as_sf()
route$Count <- ids$ccount[i]
return(route)
}
# 7. loop paths ----
get_route <- function(source,
roads = school_roads,
ids = pulled_ids){
for (i in 1:length(ids$ccount)){
route <- get_paths(i,
roads = roads,
ids = ids)
if(i == 1){
routes <- route
} else {
routes <- rbind(routes, route)
}
}
# use overline to aggregate and segmentize lines, ncores doesn't work
if(nrow(routes) == 0){
print(paste(source$nativeid, "has no routes"))
return(NULL)
} else {
routes <- routes |> overline(attrib = "Count")
# reattach nativeid to combined routes for identification
routes <- routes |> mutate(nativeid = source$nativeid)
return(routes)
}
}
# end helper functions ----
# 8. main function ----
make_school_routes <- function(id,
sources = schools,
targets = pupils,
roads = roads_clean,
WALKING_DIST = 1931.213,
BUFFER_DIST = 3218.69,
filter_counts_over = 4
){
a <- Sys.time()
# get school
school <- select_source(source_id = id,
source = sources)
school
if(nrow(school) == 0){
warning("No source found. Try a different id")
return(NULL)
}
# report id
message(paste0("Selected id: ", school$nativeid, " (", id, ")"))
message(paste0("Walking distance: ", WALKING_DIST))
message(paste0("Buffer around source for road network: ", BUFFER_DIST))
# reduce roads
message("Reducing road network...")
roads_buff <- buffer_roads(source = school,
roads = roads,
dist = BUFFER_DIST)
# get pupils for school within walking distance
school_pupils <- select_all_targets(source = school,
all_targets = targets,
dist = WALKING_DIST)
# if no pupils selected, skip to next iteration
if(nrow(school_pupils) == 0){
return(NULL)
} else {
# continue if pupils within walking distance
# blend points to roads
message("Blending points...")
school_roads <- blend_points_to_sfnetwork(source = school,
targets = school_pupils,
roads = roads_buff)
# pull node ids
message("Getting IDs...")
pulled_ids <- pull_ids(roads = school_roads)
# find routes via paths
message("Finding routes...")
school_routes <- get_route(roads = school_roads,
ids = pulled_ids,
source = school)
# filter routes by count
message(paste("Filtering routes (Count > "), filter_counts_over, ")...")
if(!is.null(school_routes)){
filtered_routes <- school_routes |>
filter(Count > filter_counts_over)
b <- Sys.time()
message(paste("Time taken:", b - a))
return(filtered_routes)
} else {
message("Returning null as no routes")
return(NULL)
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment