#' @title Isochrones calculation
#'
#' @description Calculate isochrones on a network
#'
#' @details An isochrone is the set of reachable lines around a node in a network within
#' a specified distance (or time). This function perform dynamic segmentation to return the
#' part of the edges reached and not only the fully covered edges. Several start points and
#' several distances can be given. The network can also be directed. The lines returned
#' by the function are the most accurate representation of the isochrones. However, if
#' polygons are required for mapping, the vignette "Calculating isochrones" shows
#' how to create smooth polygons from the returned sets of lines.
#'
#' @param lines A feature collection of lines representing the edges of the network
#' @param dists A vector of the size of the desired isochrones. Can also be a list of vector
#' when each start point must have its own distances. If so, the length of the list must be equal
#' to the number of rows in start_points.
#' @param start_points A feature collection of points representing the starting
#' points if the isochrones
#' @param donught A boolean indicating if the returned lines must overlap for
#' each distance (FALSE, default) or if the lines must be cut between each
#' distance step (TRUE).
#' @param mindist The minimum distance between two points. When two points are
#' too close, they might end up snapped at the same location on a line.
#' Default is 1.
#' @param weight The name of the column in lines to use an edge weight. If NULL,
#' the geographical length is used. Note that if lines are split during the
#' network creation, the weight column is recalculated proportionally to the new lines
#' length.
#' @param direction The name of the column indicating authorized
#' travelling direction on lines. if NULL, then all lines can be used in both
#' directions (undirected). The values of the column must be "FT" (From - To),
#' "TF" (To - From) or "Both".
#' @return A feature collection of lines representing the isochrones with the
#' following columns
#' \itemize{
#' \item point_id: the index of the point at the centre of the isochrone;
#' \item distance: the size of the isochrone
#' }
#' @importFrom utils strcapture
#' @export
#' @examples
#' library(sf)
#' # creating a simple network
#' wkt_lines <- c(
#' "LINESTRING (0.0 0.0, 5.0 0.0)",
#' "LINESTRING (0.0 -5.0, 5.0 -5.0)",
#' "LINESTRING (5.0 0.0, 5.0 5.0)",
#' "LINESTRING (5.0 -5.0, 5.0 -10.0)",
#' "LINESTRING (5.0 0.0, 5.0 -5.0)",
#' "LINESTRING (5.0 0.0, 10.0 0.0)",
#' "LINESTRING (5.0 -5.0, 10.0 -5.0)",
#' "LINESTRING (10.0 0, 10.0 -5.0)",
#' "LINESTRING (10.0 -10.0, 10.0 -5.0)",
#' "LINESTRING (15.0 -5.0, 10.0 -5.0)",
#' "LINESTRING (10.0 0.0, 15.0 0.0)",
#' "LINESTRING (10.0 0.0, 10.0 5.0)")
#'
#' linesdf <- data.frame(wkt = wkt_lines,
#' id = paste("l",1:length(wkt_lines),sep=""))
#'
#' lines <- st_as_sf(linesdf, wkt = "wkt", crs = 32188)
#'
#' # and the definition of the starting point
#' start_points <- data.frame(x=c(5),
#' y=c(-2.5))
#' start_points <- st_as_sf(start_points, coords = c("x","y"), crs = 32188)
#'
#' # setting the directions
#'
#' lines$direction <- "Both"
#' lines[6,"direction"] <- "TF"
#'
#' isochrones <- calc_isochrones(lines,dists = c(10,12),
#' donught = TRUE,
#' start_points = start_points,
#' direction = "direction")
calc_isochrones <- function(lines, dists, start_points, donught = FALSE, mindist = 1, weight = NULL, direction = NULL){
# if we provided a vector instead of a matrix, we create the matrix here
if(is.list(dists) == FALSE){
dist_list <- lapply(1:nrow(start_points), function(i){dists})
}else{
dist_list <- dists
}
if(length(dist_list) != nrow(start_points)){
stop("When dists is a list, it must have a length equal to the number of rows in start_points")
}
if(st_crs(start_points) != st_crs(lines)){
stop("start_points and lines must have the same crs (st_crs(start_points) == st_crs(lines))")
}
if( (is_projected(start_points) + is_projected(lines)) < 2){
stop('both lines and start_points must use a cartesian crs and not lon/lat')
}
# step1 : Check that some points are not too close to each other
# before snapping
if(nrow(start_points) >1){
xy <- st_coordinates(start_points)
#min_dists <- FNN::get.knn(xy, k = 1)
min_dists <- dbscan::kNN(xy, k = 1)
if(any(min_dists$dist < mindist)){
stop("some points are closer to each other than the mindist argument value before snapping")
}
}
# step2 : snapping the points on the lines
lines$OID <- 1:nrow(lines)
snapped_points <- snapPointsToLines2(start_points, lines, "OID")
# step3 : check the points closeness again
if(nrow(start_points) >1){
xy <- st_coordinates(snapped_points)
#min_dists <- FNN::get.knn(xy, k = 1)
min_dists <- dbscan::kNN(xy, k = 1)
if(any(min_dists$dist < mindist)){
stop("some points are closer to each other than the mindist argument value after snapping them on lines")
}
}
# step4 : splitting the lines with the points
lines$tot_length <- as.numeric(st_length(lines))
new_lines <- split_lines_at_vertex(lines, snapped_points, snapped_points$nearest_line_id, mindist = 0)
new_lines$length <- as.numeric(st_length(new_lines))
if(is.null(weight) == FALSE){
new_lines[[weight]] <- new_lines[[weight]] * (new_lines$length / new_lines$tot_length)
}else{
weight <- "length"
}
# step5 : building the graph
if (is.null(direction)){
graph_result <- build_graph(lines = new_lines, line_weight = weight, digits = 2, attrs = TRUE)
} else{
vals <- unique(new_lines[[direction]])
u <- union(vals, c("FT","TF","Both"))
if(length(u) > 3){
stop("when indicating line direction, accepted values are TF, FT and Both")
}
graph_result <- build_graph_directed(new_lines, digits = 2, line_weight = weight,
direction = direction, attrs = TRUE)
}
# finding for each start points its corresponding node in the graph
maxdist <- max(sapply(dist_list, max))
xynodes <- st_coordinates(graph_result$spvertices)
xy_points <- st_coordinates(snapped_points)
start_nodes <- dbscan::kNN(xynodes, query = xy_points, k=1)
df_start <- data.frame(
"ptOID" = 1:nrow(snapped_points),
"node_idx" = as.vector(start_nodes$id),
"node_id" = graph_result$spvertices$id[start_nodes$id],
"node_name" = graph_result$spvertices$name[start_nodes$id]
)
# creating the isochrones: sets of linestrings
all_multi_lignes <- lapply(1:nrow(df_start), function(i){
row <- df_start[i,]
dists <- dist_list[[i]]
# pour chaque point de depart calculer les distances a tous les autres points
all_dists <- t(igraph::distances(graph_result$graph,
to = row$node_idx,
v = graph_result$spvertices$id,
mode = "out"))
dist_df <- data.frame(
"node_id" = graph_result$spvertices$id,
"dist" = as.vector(all_dists)
)
# on va maintenant iterer sur toutes les distances demandees
dists2 <- c(0,dists[1:(length(dists)-1)])
all_lignes <- lapply(1:length(dists), function(j){
d <- dists[[j]]
dd <- dists2[[j]]
# on retrouve tous les noeuds a utiliser
ok_nodes <- subset(dist_df, dist_df$dist <= d)
# et toutes les edges
ok_edges <- subset(graph_result$spedges,
graph_result$spedges$start_oid %in% ok_nodes$node_id |
graph_result$spedges$end_oid %in% ok_nodes$node_id
)
# on trouve maintenant les edges sur lesquelles on doit retravailler un peu
df1 <- ok_edges
df2 <- ok_nodes
df1 <- merge(df1,df2, by.x = "start_oid",
by.y = "node_id", all.x = TRUE,
all.y = FALSE)
df1$start_dist <- df1$dist
df1$dist <- NULL
df1 <- merge(df1,df2, by.x = "end_oid",
by.y = "node_id", all.x = TRUE,
all.y = FALSE)
df1$end_dist <- df1$dist
df1$dist <- NULL
if(donught){
df1 <- subset(df1,
df1$start_dist >= dd & df1$start_dist <= d |
df1$end_dist >= dd & df1$end_dist <= d |
is.na(df1$start_dist) | is.na(df1$end_dist)
)
}
# tm_shape(lines) + tm_lines("black") + tm_shape(start_points) + tm_dots("red", size = 0.5) + tm_shape(df1) + tm_lines("blue")
# If we are in a directed graph, some edges need to be removed here
if(is.null(direction) == FALSE){
df1 <- subset(df1, (is.na(df1$end_dist) & df1$direction == "TF") == FALSE)
df1 <- subset(df1, (is.na(df1$start_dist) & df1$direction == "FT") == FALSE)
}
if(nrow(df1) == 0){
return(NULL)
}else{
ok_lines <- trim_lines_at(df1, graph_result, d, dd, i, donught)
}
# # we now have to cut the remaining edges
# test <- is.na(df1$start_dist)==FALSE & is.na(df1$end_dist)==FALSE
# no_cut <- subset(df1, test)
#
# # cutting the lines by the start
# to_cut <- subset(df1, !test)
# to_cut$node_okid <- ifelse(is.na(to_cut$start_dist), to_cut$end_oid,
# to_cut$start_oid
# )
# to_cut$ok_dist <- ifelse(is.na(to_cut$start_dist), to_cut$end_dist,
# to_cut$start_dist
# )
#
# # reordering line if required
# ext <- lines_extremities(to_cut)
# ok_nodes <- graph_result$spvertices[to_cut$node_okid,]
# test1 <- subset(ext, ext$pttype == "start")
# test2 <- subset(ext, ext$pttype == "end")
#
# d1 <- sqrt((test1$X - ok_nodes$x)**2 + (test1$Y - ok_nodes$y)**2)
# d2 <- sqrt((test2$X - ok_nodes$x)**2 + (test2$Y - ok_nodes$y)**2)
# to_keep <- subset(to_cut, d1 < d2)
# to_reverse <- subset(to_cut, d1 >= d2)
#
# # and final cut
# all_dists <- d - c(to_keep$ok_dist, to_reverse$ok_dist)
# if(nrow(to_reverse) > 0){
# all_cuts <- rbind(to_keep, reverse_lines(to_reverse))
# }else{
# all_cuts <- to_keep
# }
#
# cut_lines <- cut_lines_at_distance(all_cuts, all_dists)
#
# # saving
# ok_lines <- rbind(no_cut[c("end_oid","start_oid","edge_id","weight" )],
# cut_lines[c("end_oid","start_oid","edge_id","weight" )])
# ok_lines$distance <- d
# ok_lines$point_id <- i
# ok_lines <- ok_lines[c("point_id", "distance")]
return(ok_lines)
})
return(all_lignes[lengths(all_lignes) > 0])
})
all_multi_lignes <- do.call(rbind, unlist(all_multi_lignes, recursive = FALSE))
return(all_multi_lignes)
}
#' @title Helper for isochrones lines cutting
#'
#' @description last operation for isochrone calculation, cutting the lines at
#' their begining and ending. This is a worker function for calc_isochrones.
#'
#' @param df1 A features collection of linestrings with some specific fields.
#' @param graph_result A list produced by the functions build_graph_directed or build_graph.
#' @param d the end distance of this isochrones.
#' @param dd the start distance of this isochrones.
#' @param i the actual iteration.
#' @param donught A boolean indicating if the returned isochrone will be plained or a donught.
#'
#' @return A feature collection of lines
#'
#' @keywords internal
trim_lines_at <- function(df1, graph_result, d, dd, i, donught){
# first, we need to determine which edges will not need to be cut
if(donught){
test <- is.na(df1$start_dist) == FALSE &
is.na(df1$end_dist) == FALSE &
df1$start_dist >= dd &
df1$end_dist >= dd
}else{
test <- is.na(df1$start_dist)==FALSE & is.na(df1$end_dist)==FALSE
}
if(sum(!test) > 0){
no_cut <- subset(df1, test)
to_cut <- subset(df1, !test)
# here I need to know if the start point and end point of each line match with
# the start node and end nodes
# first, I find the coordinates of the nodes in the graph
to_cut$start_nodeX <- graph_result$spvertices[to_cut$start_oid,]$x
to_cut$end_nodeX <- graph_result$spvertices[to_cut$end_oid,]$x
to_cut$start_nodeY <- graph_result$spvertices[to_cut$start_oid,]$y
to_cut$end_nodeY <- graph_result$spvertices[to_cut$end_oid,]$y
# second, I find the coordinates of the extremities of the lines
ext <- lines_extremities(to_cut)
coords <- st_coordinates(ext)
start_coords <- subset(coords, ext$pttype == "start")
# third, I calculate the distances between the nodes and the start point
dist1 <- sqrt(((start_coords[,1] - to_cut$start_nodeX) ** 2) + ((start_coords[,2] - to_cut$start_nodeY) ** 2))
dist2 <- sqrt(((start_coords[,1] - to_cut$end_nodeX) ** 2) + ((start_coords[,2] - to_cut$end_nodeY) ** 2))
# If the distance to start_node is not the smallest, then I must reverse the geometry
# it will insure that for each line, start_node is also the first vertex of the line
test <- dist1 < dist2
to_cut2 <- rbind(
subset(to_cut, test),
reverse_lines(subset(to_cut, !test))
)
# on va passer le tout a une fonction en c++ qui va iterer et regler son cas a chaque ligne
list_of_lines <- lines_coordinates_as_list(to_cut2)
start_dists <- ifelse(is.na(to_cut2$start_dist),-1,to_cut2$start_dist)
end_dists <- ifelse(is.na(to_cut2$end_dist),-1,to_cut2$end_dist)
cut_coords <- trim_lines_for_isos_cpp(list_of_lines, start_dists, end_dists, donught, d, dd)
new_lines <- list_coordinates_as_lines(cut_coords, st_crs(to_cut2))
to_cut2$geometry <- new_lines$geometry
# saving
ok_lines <- rbind(no_cut[c("end_oid","start_oid","edge_id","weight" )],
to_cut2[c("end_oid","start_oid","edge_id","weight" )])
}else{
ok_lines <- no_cut[c("end_oid","start_oid","edge_id","weight" )]
}
ok_lines$distance <- d
ok_lines$point_id <- i
ok_lines <- ok_lines[c("point_id", "distance")]
return(ok_lines)
}
# trim_lines_at <- function(df1, graph_result, d, dd, i, donught){
#
# # we now have to cut the remaining edges
# if(donught){
# test <- is.na(df1$start_dist)==FALSE &
# is.na(df1$end_dist)==FALSE &
# df1$start_dist >= dd &
# df1$end_dist >= dd
# }else{
# test <- is.na(df1$start_dist)==FALSE & is.na(df1$end_dist)==FALSE
# }
#
# no_cut <- subset(df1, test)
# to_cut <- subset(df1, !test)
#
# if(donught){
# ok_first_cut <- is.na(to_cut$start_dist) | is.na(to_cut$end_dist)
# }else{
#
# }
#
# # cutting the lines by the start
# to_cut$node_okid <- ifelse(is.na(to_cut$start_dist), to_cut$end_oid,
# to_cut$start_oid
# )
# to_cut$ok_dist <- ifelse(is.na(to_cut$start_dist), to_cut$end_dist,
# to_cut$start_dist
# )
#
#
# # reordering line if required
# ext <- lines_extremities(to_cut)
# ok_nodes <- graph_result$spvertices[to_cut$node_okid,]
# test1 <- subset(ext, ext$pttype == "start")
# test2 <- subset(ext, ext$pttype == "end")
#
# d1 <- sqrt((test1$X - ok_nodes$x)**2 + (test1$Y - ok_nodes$y)**2)
# d2 <- sqrt((test2$X - ok_nodes$x)**2 + (test2$Y - ok_nodes$y)**2)
# to_keep <- subset(to_cut, d1 < d2)
# to_reverse <- subset(to_cut, d1 >= d2)
#
# # and final cut
# all_dists <- d - c(to_keep$ok_dist, to_reverse$ok_dist)
# if(nrow(to_reverse) > 0){
# all_cuts <- rbind(to_keep, reverse_lines(to_reverse))
# }else{
# all_cuts <- to_keep
# }
#
# cut_lines <- cut_lines_at_distance(all_cuts, all_dists)
#
# # HERE : TODO THINK ABOUT HOW TO ALSO CUT THE LINES AT THEIR BEGINING IF REQUIRED !
# if(donught){
# # HERE : TODO THINK ABOUT HOW TO ALSO CUT THE LINES AT THEIR BEGINING IF REQUIRED !
# start_cuts <- dd - cut_lines$ok_dist
# need_2nd_cut <- start_cuts > 0
# if(sum(need_2nd_cut) > 0){
# part1 <- subset(cut_lines, need_2nd_cut)
# part2 <- subset(cut_lines, !need_2nd_cut)
#
# cut_lines2 <- cut_lines_at_distance(reverse_lines(part1),
# as.vector(st_length(part1)) - start_cuts[need_2nd_cut]
# )
# cut_lines2$lineID.1 <- NULL
# cut_lines <- rbind(cut_lines2,part2)
# }
# }
#
# # saving
# ok_lines <- rbind(no_cut[c("end_oid","start_oid","edge_id","weight" )],
# cut_lines[c("end_oid","start_oid","edge_id","weight" )])
# ok_lines$distance <- d
# ok_lines$point_id <- i
# ok_lines <- ok_lines[c("point_id", "distance")]
# return(ok_lines)
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.