R/isochrones_sf.R In spNetwork: Spatial Analysis on Network

Documented in calc_isochronestrim_lines_at

#' @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)
# }


Try the spNetwork package in your browser

Any scripts or data that you put into this service are public.

spNetwork documentation built on June 22, 2024, 9:40 a.m.