R/neighbouring_sf.R
In JeremyGelb/spNetwork: Spatial Analysis on Network
Defines functions
network_listw.mc
network_listw
prepare_elements_netlistw
network_listw_worker
select_dist_function
Documented in
network_listw
network_listw.mc
network_listw_worker
prepare_elements_netlistw
select_dist_function
#' @title Select the distance to weight function
#'
#' @description Select a function to convert distance to weights
#' if a function is provided, this function will be vectorized.
#'
#' @param dist_func Could be a name in c('inverse', 'identity',
#' 'squared inverse') or a function with only one parameter x
#' @return A vectorized function used to convert distance into spatial weights
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
select_dist_function <- function(dist_func = "inverse") {
if (inherits(dist_func, "character")) {
if (dist_func == "identity") {
dist_func <- function(x) {
return(x)
}
} else if (dist_func == "inverse") {
dist_func <- function(x) {
return(1 / x)
}
} else if (dist_func == "squared inverse") {
dist_func <- function(x) {
return(1 / x^2)
}
}
}
vdist_func <- Vectorize(dist_func, "x")
return(vdist_func)
}
#defining some global variables (weird felx but ok)
utils::globalVariables(c("origin", "fid"))
# network_listw_worker_old <- function(points,lines,maxdistance,dist_func, direction=NULL, mindist=10, matrice_type = "B", verbose = FALSE, digits = 3, tol=0.1){
# ## step1 : adding the points to the lines
# lines$worker_id <- 1:nrow(lines)
# points$nearest_line_id <- as.numeric(as.character(points$nearest_line_id))
# joined <- data.table(st_drop_geometry(points))
# B <- data.table(st_drop_geometry(lines))
# joined[B, on = c("nearest_line_id" = "tmpid"),
# names(B) := mget(paste0("i.", names(B)))]
#
# if(verbose){
# print("adding the points as vertices to nearest lines")
# }
# # in the previous version, we split all the lines at their vertex
# # however, this is not required here, we could just split them
# # at points added in the network
# # new_lines <- add_vertices_lines(lines, points, joined$worker_id, 1)
# graph_lines <- split_lines_at_vertex(lines, points, joined$worker_id, 1)
#
# ## step2 : splitting the lines on vertices and adjusting weights
# if(verbose){
# print("splitting the lines for the network")
# }
# #graph_lines <- simple_lines(new_lines)
# graph_lines$lx_length <- as.numeric(sf::st_length(graph_lines))
# graph_lines$lx_weight <- (graph_lines$lx_length / graph_lines$line_length) * graph_lines$line_weight
#
# ## step3 building the network
# if(verbose){
# print("generating the network")
# }
# if (is.null(direction)){
# result_graph <- build_graph(graph_lines, digits = digits,
# attrs = TRUE, line_weight = "lx_weight")
# }else{
# #dir <- ifelse(graph_lines[[direction]]=="Both",0,1)
# result_graph <- build_graph_directed(graph_lines, digits = digits,
# attrs = TRUE, line_weight='lx_weight',
# direction = direction)
# }
# ## step4 finding for each point its corresponding vertex
# points$vertex <- closest_points(points,result_graph$spvertices)
# if(verbose){
# print("calculating the distances on the network")
# }
# starts <- subset(points,points$pttype=="start")
# u <- unique(points$vertex)
#
# ## step 4 calculating the distances between the start and end points
# base_distances <- igraph::distances(result_graph$graph,
# v = starts$vertex, to = u,
# mode = "out")
# all_ditancesdt <- data.table(base_distances)
# all_ditancesdt$origin <- starts$fid
#
# ## step 5 extracting the distances as a dataframe rather than a matrix
#
# #first groupby on rows :
# #multiple rows might correspond to the same object (if multiple points like)
# step1 <- all_ditancesdt[, lapply(.SD, min, na.rm = TRUE), by = origin ]
# originid <- step1$origin
#
# #then transpose and merge
# step2 <- transpose(step1[,2:ncol(step1)])
# step2$vertex <- u
# pts <- data.table(st_drop_geometry(points[c("fid","vertex")]))
# step2.5 <- data.table::merge.data.table(step2,pts,by="vertex",all=TRUE)
# colorder <- c(2:(ncol(step2.5)-1),1,ncol(step2.5))
# step2.5 <- data.table::setcolorder(step2.5,colorder)
#
# step3 <- step2.5[, lapply(.SD, min, na.rm=TRUE), by=fid ]
# destid <- step3$fid
# #transpose again
# step4 <- transpose(step3)
# dist_matrix <- as.matrix(step4[2:(nrow(step4)-1)])
# dist_matrix <- ifelse(is.na(dist_matrix),Inf,dist_matrix)
# rownames(dist_matrix) <- originid
# colnames(dist_matrix) <- destid
#
# ##step 6 finaly generate the neighbouring object with the data
# nblist <- lapply(1:nrow(dist_matrix),function(i){
# row <- dist_matrix[i,]
# iid <- originid[i]
# test <- row<=maxdistance & destid!=iid
# if(any(test)){
# nbs <- destid[test]
# okdistances <- row[test]
# okdistances <- ifelse(okdistances==0,mindist,okdistances)
# weights <- dist_func(okdistances)
# if(matrice_type=="B"){
# wtdweights <- rep(1,length(nbs))
# }else if (matrice_type=="W"){
# wtdweights <- weights/sum(weights)
# }else if (matrice_type == "I"){
# wtdweights <- weights
# }
# return(list(nbs,as.numeric(wtdweights)))
# }else{
# return(list(0L, NULL))
# }
# })
# nbs <- lapply(1:length(nblist), function(j) {
# nblist[[j]][[1]]
# })
# names(nbs) <- originid
# weights <- lapply(1:length(nblist), function(j) {
# nblist[[j]][[2]]
# })
# names(weights) <- originid
# return(list(nbs, weights))
# }
#' @title network_listw worker
#'
#' @description The worker function of network_listw.
#'
#' @param points A feature collection of points corresponding to start and end
#' points. It must have a column fid, grouping the points if necessary.
#' @param lines A feature collection of lines representing the network
#' @param maxdistance The maximum distance between two observation to
#' consider them as neighbours.
#' @param dist_func A vectorized function converting spatial distances into
#' weights.
#' @param mindist The minimum distance between two different observations.
#' It is important for it to be different from 0 when a W style is used.
#' @param mindist The minimum distance between two different observations.
#' It is important for it to be different from 0 when a W style is used.
#' @param direction Indicate a field giving information about authorized
#' travelling direction on lines. if NULL, then all lines can be used in both
#' directions. Must be the name of a column otherwise. The values of the
#' column must be "FT" (From - To), "TF" (To - From) or "Both".
#' @param matrice_type The type of the weighting scheme. Can be 'B' for Binary,
#' 'W' for row weighted, or 'I' (identity), see the documentation of spdep::nb2listw for details
#' @param verbose A Boolean indicating if the function should print its
#' progress
#' @param digits the number of digits to keep in the spatial coordinates (
#' simplification used to reduce risk of topological error)
#' @param tol A float indicating the spatial tolerance when points are
#' added as vertices to lines.
#' @return A list of neihbours as weights.
#' @importFrom sf st_drop_geometry
#' @importFrom utils capture.output
#' @importFrom data.table data.table .SD transpose
#' @keywords internal
#' @examples
#' #no example provided, this is an internal function
network_listw_worker<-function(points,lines,maxdistance,dist_func, direction=NULL, mindist=10, matrice_type = "B", verbose = FALSE, digits = 3, tol=0.1){
## step1 : adding the points to the lines
lines$worker_id <- 1:nrow(lines)
points$nearest_line_id <- as.numeric(as.character(points$nearest_line_id))
joined <- data.table(st_drop_geometry(points))
B <- data.table(st_drop_geometry(lines))
joined[B, on = c("nearest_line_id" = "tmpid"),
names(B) := mget(paste0("i.", names(B)))]
if(verbose){
print("adding the points as vertices to nearest lines")
}
# in the previous version, we split all the lines at their vertex
# however, this is not required here, we could just split them
# at points added in the network
# new_lines <- add_vertices_lines(lines, points, joined$worker_id, 1)
graph_lines <- split_lines_at_vertex(lines, points, joined$worker_id, mindist = mindist)
## step2 : splitting the lines on vertices and adjusting weights
if(verbose){
print("splitting the lines for the network")
}
#graph_lines <- simple_lines(new_lines)
graph_lines$lx_length <- as.numeric(sf::st_length(graph_lines))
graph_lines$lx_weight <- (graph_lines$lx_length / graph_lines$line_length) * graph_lines$line_weight
## step3 building the network
if(verbose){
print("generating the network")
}
result_graph <- build_graph_cppr(graph_lines,
digits = digits,
attrs = TRUE,
line_weight = "lx_weight",
direction = direction)
## step4 finding for each point its corresponding vertex
points$vertex <- closest_points(points,result_graph$spvertices)
points$graph_ref <- result_graph$spvertices$ref[points$vertex]
points$id_order <- 1:nrow(points)
if(verbose){
print("calculating the distances on the network")
}
starts <- subset(points,points$pttype=="start")
## step 4 calculating the distances between the start and end points
base_distances <- cppRouting::get_distance_matrix(result_graph$graph, from = starts$graph_ref, to = points$graph_ref)
base_distances[is.na(base_distances)] <- maxdistance + 1
colnames(base_distances) <- points$fid
all_ditancesdt <- data.table(base_distances)
all_ditancesdt$origin <- starts$fid
# we generate a long format for the distance matrix
all_ditancesdt_long <- data.table::melt(all_ditancesdt, id.vars = c("origin"))
# and then we can find between each feature the shortest path
ok_distances <- all_ditancesdt_long[, lapply(.SD, min, na.rm = TRUE), by = c('origin', 'variable') ]
names(ok_distances) <- c('origin', 'destination','dist')
dist_matrix <- data.table::dcast(ok_distances, origin ~ destination, value.var = 'dist')
originid <- as.integer(dist_matrix$origin)
dist_matrix$origin <- NULL
destid <- as.integer(colnames(dist_matrix))
dist_matrix <- as.matrix(dist_matrix)
##step 6 finaly generate the neighbouring object with the data
nblist <- lapply(1:nrow(dist_matrix),function(i){
row <- dist_matrix[i,]
iid <- originid[i]
test <- row<=maxdistance & destid!=iid
if(any(test)){
nbs <- destid[test]
okdistances <- row[test]
okdistances <- ifelse(okdistances==0,mindist,okdistances)
weights <- dist_func(okdistances)
if(matrice_type=="B"){
wtdweights <- rep(1,length(nbs))
}else if (matrice_type=="W"){
wtdweights <- weights/sum(weights)
}else if (matrice_type == "I"){
wtdweights <- weights
}
return(list(nbs,as.numeric(wtdweights)))
}else{
return(list(0L, NULL))
}
})
nbs <- lapply(1:length(nblist), function(j) {
nblist[[j]][[1]]
})
names(nbs) <- originid
weights <- lapply(1:length(nblist), function(j) {
nblist[[j]][[2]]
})
names(weights) <- originid
return(list(nbs, weights))
}
# prepare_elements_netlistw <- function(is,grid,snapped_points,lines,maxdistance){
#
# ## step1 preparing the spatial indices
# lines_tree <- spNetwork::build_quadtree(lines)
# snapped_points_tree <- spNetwork::build_quadtree(snapped_points)
#
# snapped_points$oids <- 1:nrow(snapped_points)
# ## step 2 extracting the quadra
# results <- lapply(is,function(i){
# quadra <- grid[i,]
# #selecting the starting points
# start_pts <- spNetwork::spatial_request(quadra,snapped_points_tree, snapped_points)
# if(nrow(start_pts)==0){
# return(NULL)
# }else{
# #start_pts <- snapped_points[snapped_points$fid %in% start_pts$fid,]
# start_pts$pttype <- "start"
# #selecting the endpoints
# poly <- st_bbox_geom(start_pts)
# buff <- st_buffer(poly, dist = maxdistance)
# if(nrow(start_pts)==nrow(snapped_points)){
# all_pts <- start_pts
# }else{
# end_pts <- spNetwork::spatial_request(buff, snapped_points_tree, snapped_points)
# end_pts <- subset(end_pts,(end_pts$oids %in% start_pts$oids)== FALSE)
# if(nrow(end_pts) > 0){
# end_pts$pttype <- "end"
# #combining all the points
# all_pts <- rbind(start_pts,end_pts)
# }else{
# all_pts <- start_pts
# }
#
# }
# #selecting the lines
# selected_lines <- spNetwork::spatial_request(buff, lines_tree, lines)
# #calculating the elements
# return(list(all_pts,selected_lines))
# }
# })
# return(results)
# }
#' @title Data preparation for network_listw
#'
#' @description Function to prepare selected points and selected lines during
#' the process.
#'
#' @param is The indices of the quadras to use in the grid
#' @param grid A feature collection of polygons representing the quadras to split
#' calculus
#' @param snapped_points The start and end points snapped to the lines
#' @param lines The lines representing the network
#' @param maxdistance The maximum distance between two observation to considere
#' them as neighbours.
#' @return A list of two elements : selected points and selected lines
#' @keywords internal
#' @importFrom sf st_bbox st_buffer
#' @examples
#' #no example provided, this is an internal function
prepare_elements_netlistw <- function(is,grid,snapped_points,lines,maxdistance){
## step1 preparing the spatial indices
grid$grid_id <- as.character(is)
grid$oid <- NULL
# We precalculate the intersections
inter_pts <- st_join(snapped_points, grid)
inter_pts <- split(inter_pts, inter_pts$grid_id)
inter_end_pts <- st_join(snapped_points, st_buffer(grid, maxdistance))
inter_end_pts <- split(inter_end_pts, inter_end_pts$grid_id)
inter_lines <- st_join(lines, st_buffer(grid, maxdistance))
inter_lines <- split(inter_lines, inter_lines$grid_id)
snapped_points$oids <- 1:nrow(snapped_points)
## step 2 extracting the quadra
results <- lapply(is,function(i){
quadra <- grid[i,]
gid <- quadra$grid_id
#selecting the starting points
start_pts <- inter_pts[[gid]]
if(is.null(start_pts)){
return(NULL)
}
if(nrow(start_pts)==0){
return(NULL)
}else{
#start_pts <- snapped_points[snapped_points$fid %in% start_pts$fid,]
start_pts$pttype <- "start"
#selecting the endpoints
if(nrow(start_pts)==nrow(snapped_points)){
all_pts <- start_pts
}else{
end_pts <- inter_end_pts[[gid]]
end_pts <- subset(end_pts,(end_pts$oids %in% start_pts$oids)== FALSE)
if(nrow(end_pts) > 0){
end_pts$pttype <- "end"
#combining all the points
all_pts <- rbind(start_pts,end_pts)
}else{
all_pts <- start_pts
}
}
#selecting the lines
selected_lines <- inter_lines[[gid]]
#calculating the elements
return(list(all_pts,selected_lines))
}
})
return(results)
}
#' @title Network distance listw
#'
#' @description Generate listw object (spdep like) based on network distances.
#'
#' @param origins A feature collection of lines, points, or polygons
#' for which the spatial neighbouring list will be built
#' @param lines A feature collection of lines representing the network
#' @param maxdistance The maximum distance between two observations to
#' consider them as neighbours.
#' @param method A string indicating how the starting points will be built.
#' If 'centroid' is used, then the centre of lines or polygons is used. If
#' 'pointsalong' is used, then points will be placed along polygons' borders or
#' along lines as starting and end points. If 'ends' is used (only for lines)
#' the first and last vertices of lines are used as starting and ending points.
#' @param point_dist A float, defining the distance between points when the
#' method 'pointsalong' is selected.
#' @param snap_dist The maximum distance to snap the start and end points on
#' the network.
#' @param line_weight The weighting to use for lines. Default is "length"
#' (the geographical length), but can be the name of a column. The value is
#' considered proportional to the geographical length of the lines.
#' @param mindist The minimum distance between two different observations.
#' It is important for it to be different from 0 when a W style is used.
#' @param direction Indicates a field providing information about authorized
#' travelling direction on lines. if NULL, then all lines can be used in both
#' directions. Must be the name of a column otherwise. The values of the
#' column must be "FT" (From - To), "TF" (To - From) or "Both".
#' @param dist_func Indicates the function to use to convert the distance
#' between observation in spatial weights. Can be 'identity', 'inverse',
#' 'squared inverse' or a function with one parameter x that will be
#' vectorized internally
#' @param matrice_type The type of the weighting scheme. Can be 'B' for Binary,
#' 'W' for row weighted, or 'I' (identity), see the documentation of spdep::nb2listw for details
#' @param grid_shape A vector of length 2 indicating the shape of the grid to
#' use for splitting the dataset. Default is c(1,1), so all the calculation is
#' done in one go. It might be necessary to split it if the dataset is large.
#' @param verbose A Boolean indicating if the function should print its
#' progress
#' @param digits The number of digits to retain in the spatial coordinates (
#' simplification used to reduce risk of topological error)
#' @param tol A float indicating the spatial tolerance when points are
#' added as vertices to lines.
#' @return A listw object (spdep like) if matrice_type is "B" or "W". If
#' matrice_type is I, then a list with a nblist object and a list of weights
#' is returned.
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom graphics plot
#' @importFrom sf st_point_on_surface st_length
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' listw <- network_listw(mtl_network,
#' mtl_network,
#' maxdistance = 500,
#' method = "centroid",
#' line_weight = "length",
#' dist_func = 'squared inverse',
#' matrice_type='B',
#' grid_shape = c(2,2))
#' }
network_listw <- function(origins,lines, maxdistance, method="centroid", point_dist=NULL, snap_dist=Inf, line_weight = "length", mindist=10, direction=NULL, dist_func = "inverse", matrice_type = "B", grid_shape=c(1,1), verbose = FALSE, digits = 3, tol=0.1){
## step1 adjusting the weights of the lines
lines$line_length <- as.numeric(sf::st_length(lines))
if(line_weight=="length"){
lines$line_weight <- as.numeric(sf::st_length(lines))
}else {
lines$line_weight <- lines[[line_weight]]
}
if(min(lines$line_weight)<=0){
stop("the weights of the lines must be superior to 0")
}
## step2 adjusting the directions of the lines (done now by the graph function)
#if(is.null(direction) == FALSE){
# lines <- lines_direction(lines,direction)
#}
## step3 checking the matrix type
if (matrice_type %in% c("B", "W", "I") == FALSE) {
stop("Matrice type must be B (binary), W (row standardized) or I (identity)")
}
## step 4 creating the vectorized distance function
vdist_func <- select_dist_function(dist_func)
lines$tmpid <- 1:nrow(lines)
## step5 generating the starting points
if(verbose){
print("generating the starting points")
}
origins$fid <- 1:nrow(origins)
if(unique(st_geometry_type(origins))=="POLYGON"){
if(method=="centroid"){
centers_geom <- st_point_on_surface(origins$geometry)
centers <- origins
centers$geometry <- centers_geom
}else if(method=="pointsalong"){
centers <- surrounding_points(origins,point_dist)
}
}else if(unique(st_geometry_type(origins))=="POINT"){
centers <- origins
}else if(unique(st_geometry_type(origins))=="LINESTRING"){
if(method=="centroid"){
if(verbose){
print("getting the centers of the lines ...")
centers <- lines_center(origins)
}else{
invisible(capture.output(centers <- lines_center(origins)))
}
}else if(method=="ends"){
centers <- lines_extremities(origins)
}else if(method=="pointsalong"){
centers <- lines_points_along(origins,point_dist)
}else{
stop("if origins class is feature collection of linestrings, method must be on of centroid, ends or pointsalong")
}
}else{
stop("origins must be a feature collection of points, polygons or linestrings")
}
if(verbose){
print("snapping the points to the lines (only once)")
}
#snapped_points <- maptools::snapPointsToLines(centers,lines,maxDist = snap_dist, idField="tmpid")
snapped_points <- snapPointsToLines2(centers,lines, idField="tmpid")
#snapped_points <- cbind(snapped_points, centers)
## step 6 building grid
mygrid <- build_grid(grid_shape,list(origins,lines))
if(verbose){
print("starting the network part")
}
ids <- 1:nrow(mygrid)
list_elements <- prepare_elements_netlistw(ids,mygrid,snapped_points,lines,maxdistance)
## step7 iterating over the grid
listvalues <- lapply(1:nrow(mygrid),function(i){
if(verbose){
print(paste("working on quadra : ",i,"/",max(ids),sep=""))
}
elements <- list_elements[[i]]
if(length(elements)==0){
return()
}else {
all_pts <- elements[[1]]
selected_lines <- elements[[2]]
#calculating the elements
values <- network_listw_worker(points = all_pts,
lines = selected_lines,
maxdistance = maxdistance,
direction=direction,
mindist=mindist,
dist_func = vdist_func,
matrice_type = matrice_type,
verbose = verbose,
digits = digits,
tol=tol)
return(values)
}
})
if(verbose){
print("building the final listw object")
}
## step8 combining the results
okvalues <- listvalues[lengths(listvalues) != 0]
nblist <- do.call("c", lapply(okvalues, function(value) {
return(value[[1]])
}))
ordered_nblist <- lapply(as.character(origins$fid), function(x) {
return(nblist[[x]])
})
nbweights <- do.call("c", lapply(okvalues, function(value) {
return(value[[2]])
}))
names(nbweights) <- names(nblist)
ordered_weights <- lapply(as.character(origins$fid), function(x) {
return(nbweights[[x]])
})
# setting the nb attributes
attr(ordered_nblist, "class") <- "nb"
attr(ordered_nblist, "region.id") <- as.character(origins$fid)
if(verbose){
print("finally generating the listw object ...")
}
if (matrice_type != "I"){
# setting the final listw attributes
listw <- spdep::nb2listw(ordered_nblist, glist = ordered_weights,
zero.policy = TRUE, style = matrice_type)
return(listw)
}else{
return(list("nb_list" = ordered_nblist,
"weights" = ordered_weights))
}
}
#' @title Network distance listw (multicore)
#'
#' @description Generate listw object (spdep like) based on network distances with multicore support.
#'
#' @param origins A feature collection of linestrings, points or polygons for
#' which the spatial neighbouring list will be built.
#' @param lines A feature collection of linestrings representing the network
#' @param maxdistance The maximum distance between two observations to consider
#' them as neighbours.
#' @param method A string indicating how the starting points will be built. If
#' 'centroid' is used, then the centre of lines or polygons is used. If
#' 'pointsalong' is used, then points will be placed along polygons' borders
#' or along lines as starting and end points. If 'ends' is used (only for
#' lines) the first and last vertices of lines are used as starting and ending
#' points.
#' @param point_dist A float, defining the distance between points when the
#' method pointsalong is selected.
#' @param snap_dist the maximum distance to snap the start and end points on the
#' network.
#' @param line_weight The weights to use for lines. Default is "length" (the
#' geographical length), but can be the name of a column. The value is
#' considered proportional with the geographical length of the lines.
#' @param mindist The minimum distance between two different observations. It is
#' important for it to be different from 0 when a W style is used.
#' @param direction Indicates a field giving information about authorized
#' travelling direction on lines. if NULL, then all lines can be used in both
#' directions. Must be the name of a column otherwise. The values of the
#' column must be "FT" (From - To), "TF" (To - From) or "Both".
#' @param dist_func Indicates the function to use to convert the distance
#' between observation in spatial weights. Can be 'identity', 'inverse',
#' 'squared inverse' or a function with one parameter x that will be
#' vectorized internally
#' @param matrice_type The type of the weighting scheme. Can be 'B' for Binary,
#' 'W' for row weighted, or 'I' (identity) see the documentation of
#' spdep::nb2listw for details
#' @param grid_shape A vector of length 2 indicating the shape of the grid to
#' use for splitting the dataset. Default is c(1,1), so all the calculation is
#' done in one go. It might be necessary to split it if the dataset is large.
#' @param verbose A Boolean indicating if the function should print its progress
#' @param digits The number of digits to retain in the spatial coordinates (
#' simplification used to reduce risk of topological error)
#' @param tol A float indicating the spatial tolerance when points are added as
#' vertices to lines.
#' @return A listw object (spdep like) if matrice_type is "B" or "W". If
#' matrice_type is I, then a list with a nblist object and a list of weights
#' is returned.
#' @importFrom sf st_length st_point_on_surface
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' future::plan(future::multisession(workers=1))
#' listw <- network_listw.mc(mtl_network,mtl_network,maxdistance=500,
#' method = "centroid", line_weight = "length",
#' dist_func = 'squared inverse', matrice_type='B', grid_shape = c(2,2))
#' ## make sure any open connections are closed afterward
#' if (!inherits(future::plan(), "sequential")) future::plan(future::sequential)
#'}
network_listw.mc <- function(origins,lines,maxdistance, method="centroid", point_dist=NULL, snap_dist=Inf, line_weight = "length", mindist=10, direction=NULL, dist_func = "inverse", matrice_type = "B", grid_shape=c(1,1), verbose = FALSE, digits = 3, tol=0.1){
##adjusting the weights of the lines
lines$line_length <- as.numeric(st_length(lines))
if(line_weight=="length"){
lines$line_weight <- as.numeric(st_length(lines))
}else {
lines$line_weight <- lines[[line_weight]]
}
if(min(lines$line_weight)<=0){
stop("the weights of the lines must be superior to 0")
}
# adjusting the directions of the lines
if(is.null(direction) == FALSE){
lines <- lines_direction(lines,direction)
}
## checking the matrix type
if (matrice_type %in% c("B", "W", "I") == FALSE) {
stop("Matrice type must be B (binary), W (row standardized) or I (identity)")
}
## creating the vectorized distance function
vdist_func <- select_dist_function(dist_func)
## setting the OID : tmpid
lines$tmpid <- 1:nrow(lines)
##generating the starting points
## step5 generating the starting points
if(verbose){
print("generating the starting points")
}
origins$fid <- 1:nrow(origins)
if(unique(st_geometry_type(origins))=="POLYGON"){
if(method=="centroid"){
centers_geom <- st_point_on_surface(origins$geometry)
centers <- origins
centers$geometry <- centers_geom
}else if(method=="pointsalong"){
centers <- surrounding_points(origins,point_dist)
}
}else if(unique(st_geometry_type(origins))=="POINT"){
centers <- origins
}else if(unique(st_geometry_type(origins))=="LINESTRING"){
if(method=="centroid"){
if(verbose){
print("getting the centers of the lines ...")
centers <- lines_center(origins)
}else{
invisible(capture.output(centers <- lines_center(origins)))
}
}else if(method=="ends"){
centers <- lines_extremities(origins)
}else if(method=="pointsalong"){
centers <- lines_points_along(origins,point_dist)
}else{
stop("if origins class is feature collection of linestrings, method must be on of centroid, ends or pointsalong")
}
}else{
stop("origins must be a feature collection of points, polygons or linestrings")
}
if(verbose){
print("snapping the points to the lines (only once)")
}
#snapped_points <- maptools::snapPointsToLines(centers,lines,maxDist = snap_dist, idField="tmpid")
snapped_points <- snapPointsToLines2(centers,lines, idField="tmpid")
#snapped_points <- cbind(snapped_points, centers)
##building grid
grid <- build_grid(grid_shape,list(origins,lines))
if(verbose){
print("starting the network part")
}
all_is <- 1:nrow(grid)
iseq <- list()
cnt <- 0
for(i in 1:grid_shape[[1]]){
start <- cnt*grid_shape[[2]]+1
iseq[[length(iseq)+1]] <- list(cnt+1,all_is[start:(start+grid_shape[[2]]-1)])
cnt<-cnt+1
}
listelements <- future.apply::future_lapply(iseq,function(ii){
elements <- prepare_elements_netlistw(ii[[2]],grid,snapped_points,lines,maxdistance)
return(elements)
}, future.packages = c("sf", "spNetwork"))
listelements <- unlist(listelements,recursive = FALSE)
##iterating on the grid
listvalues <- future.apply::future_lapply(listelements,function(elements){
##step1 : preparing elements
if(is.null(elements)){
return()
}else {
all_pts <- elements[[1]]
selected_lines <- elements[[2]]
#calculating the elements
values <- network_listw_worker(all_pts, selected_lines,
maxdistance, direction=direction, mindist=mindist,
dist_func = vdist_func, matrice_type = matrice_type,
verbose = verbose, digits = digits, tol=tol)
return(values)
}
}, future.packages = c("sf", "spNetwork"))
if(verbose){
print("building the final listw object")
}
okvalues <- listvalues[lengths(listvalues) != 0]
nblist <- do.call("c", lapply(okvalues, function(value) {
return(value[[1]])
}))
ordered_nblist <- lapply(as.character(origins$fid), function(x) {
return(nblist[[x]])
})
nbweights <- do.call("c", lapply(okvalues, function(value) {
return(value[[2]])
}))
names(nbweights) <- names(nblist)
ordered_weights <- lapply(as.character(origins$fid), function(x) {
return(nbweights[[x]])
})
## setting the nb attributes
attr(ordered_nblist, "class") <- "nb"
attr(ordered_nblist, "region.id") <- as.character(origins$fid)
if(verbose){
print("finally generating the listw object ...")
}
if (matrice_type != "I"){
# setting the final listw attributes
listw <- spdep::nb2listw(ordered_nblist, glist = ordered_weights,
zero.policy = TRUE, style = matrice_type)
return(listw)
}else{
return(list("nb_list" = ordered_nblist,
"weights" = ordered_weights))
}
}
JeremyGelb/spNetwork documentation built on July 4, 2025, 1:50 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions network_listw.mc network_listw prepare_elements_netlistw network_listw_worker select_dist_function
Documented in network_listw network_listw.mc network_listw_worker prepare_elements_netlistw select_dist_function
#' @title Select the distance to weight function
#'
#' @description Select a function to convert distance to weights
#' if a function is provided, this function will be vectorized.
#'
#' @param dist_func Could be a name in c('inverse', 'identity',
#' 'squared inverse') or a function with only one parameter x
#' @return A vectorized function used to convert distance into spatial weights
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
select_dist_function <- function(dist_func = "inverse") {
if (inherits(dist_func, "character")) {
if (dist_func == "identity") {
dist_func <- function(x) {
return(x)
}
} else if (dist_func == "inverse") {
dist_func <- function(x) {
return(1 / x)
}
} else if (dist_func == "squared inverse") {
dist_func <- function(x) {
return(1 / x^2)
}
}
}
vdist_func <- Vectorize(dist_func, "x")
return(vdist_func)
}
#defining some global variables (weird felx but ok)
utils::globalVariables(c("origin", "fid"))
# network_listw_worker_old <- function(points,lines,maxdistance,dist_func, direction=NULL, mindist=10, matrice_type = "B", verbose = FALSE, digits = 3, tol=0.1){
# ## step1 : adding the points to the lines
# lines$worker_id <- 1:nrow(lines)
# points$nearest_line_id <- as.numeric(as.character(points$nearest_line_id))
# joined <- data.table(st_drop_geometry(points))
# B <- data.table(st_drop_geometry(lines))
# joined[B, on = c("nearest_line_id" = "tmpid"),
# names(B) := mget(paste0("i.", names(B)))]
#
# if(verbose){
# print("adding the points as vertices to nearest lines")
# }
# # in the previous version, we split all the lines at their vertex
# # however, this is not required here, we could just split them
# # at points added in the network
# # new_lines <- add_vertices_lines(lines, points, joined$worker_id, 1)
# graph_lines <- split_lines_at_vertex(lines, points, joined$worker_id, 1)
#
# ## step2 : splitting the lines on vertices and adjusting weights
# if(verbose){
# print("splitting the lines for the network")
# }
# #graph_lines <- simple_lines(new_lines)
# graph_lines$lx_length <- as.numeric(sf::st_length(graph_lines))
# graph_lines$lx_weight <- (graph_lines$lx_length / graph_lines$line_length) * graph_lines$line_weight
#
# ## step3 building the network
# if(verbose){
# print("generating the network")
# }
# if (is.null(direction)){
# result_graph <- build_graph(graph_lines, digits = digits,
# attrs = TRUE, line_weight = "lx_weight")
# }else{
# #dir <- ifelse(graph_lines[[direction]]=="Both",0,1)
# result_graph <- build_graph_directed(graph_lines, digits = digits,
# attrs = TRUE, line_weight='lx_weight',
# direction = direction)
# }
# ## step4 finding for each point its corresponding vertex
# points$vertex <- closest_points(points,result_graph$spvertices)
# if(verbose){
# print("calculating the distances on the network")
# }
# starts <- subset(points,points$pttype=="start")
# u <- unique(points$vertex)
#
# ## step 4 calculating the distances between the start and end points
# base_distances <- igraph::distances(result_graph$graph,
# v = starts$vertex, to = u,
# mode = "out")
# all_ditancesdt <- data.table(base_distances)
# all_ditancesdt$origin <- starts$fid
#
# ## step 5 extracting the distances as a dataframe rather than a matrix
#
# #first groupby on rows :
# #multiple rows might correspond to the same object (if multiple points like)
# step1 <- all_ditancesdt[, lapply(.SD, min, na.rm = TRUE), by = origin ]
# originid <- step1$origin
#
# #then transpose and merge
# step2 <- transpose(step1[,2:ncol(step1)])
# step2$vertex <- u
# pts <- data.table(st_drop_geometry(points[c("fid","vertex")]))
# step2.5 <- data.table::merge.data.table(step2,pts,by="vertex",all=TRUE)
# colorder <- c(2:(ncol(step2.5)-1),1,ncol(step2.5))
# step2.5 <- data.table::setcolorder(step2.5,colorder)
#
# step3 <- step2.5[, lapply(.SD, min, na.rm=TRUE), by=fid ]
# destid <- step3$fid
# #transpose again
# step4 <- transpose(step3)
# dist_matrix <- as.matrix(step4[2:(nrow(step4)-1)])
# dist_matrix <- ifelse(is.na(dist_matrix),Inf,dist_matrix)
# rownames(dist_matrix) <- originid
# colnames(dist_matrix) <- destid
#
# ##step 6 finaly generate the neighbouring object with the data
# nblist <- lapply(1:nrow(dist_matrix),function(i){
# row <- dist_matrix[i,]
# iid <- originid[i]
# test <- row<=maxdistance & destid!=iid
# if(any(test)){
# nbs <- destid[test]
# okdistances <- row[test]
# okdistances <- ifelse(okdistances==0,mindist,okdistances)
# weights <- dist_func(okdistances)
# if(matrice_type=="B"){
# wtdweights <- rep(1,length(nbs))
# }else if (matrice_type=="W"){
# wtdweights <- weights/sum(weights)
# }else if (matrice_type == "I"){
# wtdweights <- weights
# }
# return(list(nbs,as.numeric(wtdweights)))
# }else{
# return(list(0L, NULL))
# }
# })
# nbs <- lapply(1:length(nblist), function(j) {
# nblist[[j]][[1]]
# })
# names(nbs) <- originid
# weights <- lapply(1:length(nblist), function(j) {
# nblist[[j]][[2]]
# })
# names(weights) <- originid
# return(list(nbs, weights))
# }
#' @title network_listw worker
#'
#' @description The worker function of network_listw.
#'
#' @param points A feature collection of points corresponding to start and end
#' points. It must have a column fid, grouping the points if necessary.
#' @param lines A feature collection of lines representing the network
#' @param maxdistance The maximum distance between two observation to
#' consider them as neighbours.
#' @param dist_func A vectorized function converting spatial distances into
#' weights.
#' @param mindist The minimum distance between two different observations.
#' It is important for it to be different from 0 when a W style is used.
#' @param mindist The minimum distance between two different observations.
#' It is important for it to be different from 0 when a W style is used.
#' @param direction Indicate a field giving information about authorized
#' travelling direction on lines. if NULL, then all lines can be used in both
#' directions. Must be the name of a column otherwise. The values of the
#' column must be "FT" (From - To), "TF" (To - From) or "Both".
#' @param matrice_type The type of the weighting scheme. Can be 'B' for Binary,
#' 'W' for row weighted, or 'I' (identity), see the documentation of spdep::nb2listw for details
#' @param verbose A Boolean indicating if the function should print its
#' progress
#' @param digits the number of digits to keep in the spatial coordinates (
#' simplification used to reduce risk of topological error)
#' @param tol A float indicating the spatial tolerance when points are
#' added as vertices to lines.
#' @return A list of neihbours as weights.
#' @importFrom sf st_drop_geometry
#' @importFrom utils capture.output
#' @importFrom data.table data.table .SD transpose
#' @keywords internal
#' @examples
#' #no example provided, this is an internal function
network_listw_worker<-function(points,lines,maxdistance,dist_func, direction=NULL, mindist=10, matrice_type = "B", verbose = FALSE, digits = 3, tol=0.1){
## step1 : adding the points to the lines
lines$worker_id <- 1:nrow(lines)
points$nearest_line_id <- as.numeric(as.character(points$nearest_line_id))
joined <- data.table(st_drop_geometry(points))
B <- data.table(st_drop_geometry(lines))
joined[B, on = c("nearest_line_id" = "tmpid"),
names(B) := mget(paste0("i.", names(B)))]
if(verbose){
print("adding the points as vertices to nearest lines")
}
# in the previous version, we split all the lines at their vertex
# however, this is not required here, we could just split them
# at points added in the network
# new_lines <- add_vertices_lines(lines, points, joined$worker_id, 1)
graph_lines <- split_lines_at_vertex(lines, points, joined$worker_id, mindist = mindist)
## step2 : splitting the lines on vertices and adjusting weights
if(verbose){
print("splitting the lines for the network")
}
#graph_lines <- simple_lines(new_lines)
graph_lines$lx_length <- as.numeric(sf::st_length(graph_lines))
graph_lines$lx_weight <- (graph_lines$lx_length / graph_lines$line_length) * graph_lines$line_weight
## step3 building the network
if(verbose){
print("generating the network")
}
result_graph <- build_graph_cppr(graph_lines,
digits = digits,
attrs = TRUE,
line_weight = "lx_weight",
direction = direction)
## step4 finding for each point its corresponding vertex
points$vertex <- closest_points(points,result_graph$spvertices)
points$graph_ref <- result_graph$spvertices$ref[points$vertex]
points$id_order <- 1:nrow(points)
if(verbose){
print("calculating the distances on the network")
}
starts <- subset(points,points$pttype=="start")
## step 4 calculating the distances between the start and end points
base_distances <- cppRouting::get_distance_matrix(result_graph$graph, from = starts$graph_ref, to = points$graph_ref)
base_distances[is.na(base_distances)] <- maxdistance + 1
colnames(base_distances) <- points$fid
all_ditancesdt <- data.table(base_distances)
all_ditancesdt$origin <- starts$fid
# we generate a long format for the distance matrix
all_ditancesdt_long <- data.table::melt(all_ditancesdt, id.vars = c("origin"))
# and then we can find between each feature the shortest path
ok_distances <- all_ditancesdt_long[, lapply(.SD, min, na.rm = TRUE), by = c('origin', 'variable') ]
names(ok_distances) <- c('origin', 'destination','dist')
dist_matrix <- data.table::dcast(ok_distances, origin ~ destination, value.var = 'dist')
originid <- as.integer(dist_matrix$origin)
dist_matrix$origin <- NULL
destid <- as.integer(colnames(dist_matrix))
dist_matrix <- as.matrix(dist_matrix)
##step 6 finaly generate the neighbouring object with the data
nblist <- lapply(1:nrow(dist_matrix),function(i){
row <- dist_matrix[i,]
iid <- originid[i]
test <- row<=maxdistance & destid!=iid
if(any(test)){
nbs <- destid[test]
okdistances <- row[test]
okdistances <- ifelse(okdistances==0,mindist,okdistances)
weights <- dist_func(okdistances)
if(matrice_type=="B"){
wtdweights <- rep(1,length(nbs))
}else if (matrice_type=="W"){
wtdweights <- weights/sum(weights)
}else if (matrice_type == "I"){
wtdweights <- weights
}
return(list(nbs,as.numeric(wtdweights)))
}else{
return(list(0L, NULL))
}
})
nbs <- lapply(1:length(nblist), function(j) {
nblist[[j]][[1]]
})
names(nbs) <- originid
weights <- lapply(1:length(nblist), function(j) {
nblist[[j]][[2]]
})
names(weights) <- originid
return(list(nbs, weights))
}
# prepare_elements_netlistw <- function(is,grid,snapped_points,lines,maxdistance){
#
# ## step1 preparing the spatial indices
# lines_tree <- spNetwork::build_quadtree(lines)
# snapped_points_tree <- spNetwork::build_quadtree(snapped_points)
#
# snapped_points$oids <- 1:nrow(snapped_points)
# ## step 2 extracting the quadra
# results <- lapply(is,function(i){
# quadra <- grid[i,]
# #selecting the starting points
# start_pts <- spNetwork::spatial_request(quadra,snapped_points_tree, snapped_points)
# if(nrow(start_pts)==0){
# return(NULL)
# }else{
# #start_pts <- snapped_points[snapped_points$fid %in% start_pts$fid,]
# start_pts$pttype <- "start"
# #selecting the endpoints
# poly <- st_bbox_geom(start_pts)
# buff <- st_buffer(poly, dist = maxdistance)
# if(nrow(start_pts)==nrow(snapped_points)){
# all_pts <- start_pts
# }else{
# end_pts <- spNetwork::spatial_request(buff, snapped_points_tree, snapped_points)
# end_pts <- subset(end_pts,(end_pts$oids %in% start_pts$oids)== FALSE)
# if(nrow(end_pts) > 0){
# end_pts$pttype <- "end"
# #combining all the points
# all_pts <- rbind(start_pts,end_pts)
# }else{
# all_pts <- start_pts
# }
#
# }
# #selecting the lines
# selected_lines <- spNetwork::spatial_request(buff, lines_tree, lines)
# #calculating the elements
# return(list(all_pts,selected_lines))
# }
# })
# return(results)
# }
#' @title Data preparation for network_listw
#'
#' @description Function to prepare selected points and selected lines during
#' the process.
#'
#' @param is The indices of the quadras to use in the grid
#' @param grid A feature collection of polygons representing the quadras to split
#' calculus
#' @param snapped_points The start and end points snapped to the lines
#' @param lines The lines representing the network
#' @param maxdistance The maximum distance between two observation to considere
#' them as neighbours.
#' @return A list of two elements : selected points and selected lines
#' @keywords internal
#' @importFrom sf st_bbox st_buffer
#' @examples
#' #no example provided, this is an internal function
prepare_elements_netlistw <- function(is,grid,snapped_points,lines,maxdistance){
## step1 preparing the spatial indices
grid$grid_id <- as.character(is)
grid$oid <- NULL
# We precalculate the intersections
inter_pts <- st_join(snapped_points, grid)
inter_pts <- split(inter_pts, inter_pts$grid_id)
inter_end_pts <- st_join(snapped_points, st_buffer(grid, maxdistance))
inter_end_pts <- split(inter_end_pts, inter_end_pts$grid_id)
inter_lines <- st_join(lines, st_buffer(grid, maxdistance))
inter_lines <- split(inter_lines, inter_lines$grid_id)
snapped_points$oids <- 1:nrow(snapped_points)
## step 2 extracting the quadra
results <- lapply(is,function(i){
quadra <- grid[i,]
gid <- quadra$grid_id
#selecting the starting points
start_pts <- inter_pts[[gid]]
if(is.null(start_pts)){
return(NULL)
}
if(nrow(start_pts)==0){
return(NULL)
}else{
#start_pts <- snapped_points[snapped_points$fid %in% start_pts$fid,]
start_pts$pttype <- "start"
#selecting the endpoints
if(nrow(start_pts)==nrow(snapped_points)){
all_pts <- start_pts
}else{
end_pts <- inter_end_pts[[gid]]
end_pts <- subset(end_pts,(end_pts$oids %in% start_pts$oids)== FALSE)
if(nrow(end_pts) > 0){
end_pts$pttype <- "end"
#combining all the points
all_pts <- rbind(start_pts,end_pts)
}else{
all_pts <- start_pts
}
}
#selecting the lines
selected_lines <- inter_lines[[gid]]
#calculating the elements
return(list(all_pts,selected_lines))
}
})
return(results)
}
#' @title Network distance listw
#'
#' @description Generate listw object (spdep like) based on network distances.
#'
#' @param origins A feature collection of lines, points, or polygons
#' for which the spatial neighbouring list will be built
#' @param lines A feature collection of lines representing the network
#' @param maxdistance The maximum distance between two observations to
#' consider them as neighbours.
#' @param method A string indicating how the starting points will be built.
#' If 'centroid' is used, then the centre of lines or polygons is used. If
#' 'pointsalong' is used, then points will be placed along polygons' borders or
#' along lines as starting and end points. If 'ends' is used (only for lines)
#' the first and last vertices of lines are used as starting and ending points.
#' @param point_dist A float, defining the distance between points when the
#' method 'pointsalong' is selected.
#' @param snap_dist The maximum distance to snap the start and end points on
#' the network.
#' @param line_weight The weighting to use for lines. Default is "length"
#' (the geographical length), but can be the name of a column. The value is
#' considered proportional to the geographical length of the lines.
#' @param mindist The minimum distance between two different observations.
#' It is important for it to be different from 0 when a W style is used.
#' @param direction Indicates a field providing information about authorized
#' travelling direction on lines. if NULL, then all lines can be used in both
#' directions. Must be the name of a column otherwise. The values of the
#' column must be "FT" (From - To), "TF" (To - From) or "Both".
#' @param dist_func Indicates the function to use to convert the distance
#' between observation in spatial weights. Can be 'identity', 'inverse',
#' 'squared inverse' or a function with one parameter x that will be
#' vectorized internally
#' @param matrice_type The type of the weighting scheme. Can be 'B' for Binary,
#' 'W' for row weighted, or 'I' (identity), see the documentation of spdep::nb2listw for details
#' @param grid_shape A vector of length 2 indicating the shape of the grid to
#' use for splitting the dataset. Default is c(1,1), so all the calculation is
#' done in one go. It might be necessary to split it if the dataset is large.
#' @param verbose A Boolean indicating if the function should print its
#' progress
#' @param digits The number of digits to retain in the spatial coordinates (
#' simplification used to reduce risk of topological error)
#' @param tol A float indicating the spatial tolerance when points are
#' added as vertices to lines.
#' @return A listw object (spdep like) if matrice_type is "B" or "W". If
#' matrice_type is I, then a list with a nblist object and a list of weights
#' is returned.
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom graphics plot
#' @importFrom sf st_point_on_surface st_length
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' listw <- network_listw(mtl_network,
#' mtl_network,
#' maxdistance = 500,
#' method = "centroid",
#' line_weight = "length",
#' dist_func = 'squared inverse',
#' matrice_type='B',
#' grid_shape = c(2,2))
#' }
network_listw <- function(origins,lines, maxdistance, method="centroid", point_dist=NULL, snap_dist=Inf, line_weight = "length", mindist=10, direction=NULL, dist_func = "inverse", matrice_type = "B", grid_shape=c(1,1), verbose = FALSE, digits = 3, tol=0.1){
## step1 adjusting the weights of the lines
lines$line_length <- as.numeric(sf::st_length(lines))
if(line_weight=="length"){
lines$line_weight <- as.numeric(sf::st_length(lines))
}else {
lines$line_weight <- lines[[line_weight]]
}
if(min(lines$line_weight)<=0){
stop("the weights of the lines must be superior to 0")
}
## step2 adjusting the directions of the lines (done now by the graph function)
#if(is.null(direction) == FALSE){
# lines <- lines_direction(lines,direction)
#}
## step3 checking the matrix type
if (matrice_type %in% c("B", "W", "I") == FALSE) {
stop("Matrice type must be B (binary), W (row standardized) or I (identity)")
}
## step 4 creating the vectorized distance function
vdist_func <- select_dist_function(dist_func)
lines$tmpid <- 1:nrow(lines)
## step5 generating the starting points
if(verbose){
print("generating the starting points")
}
origins$fid <- 1:nrow(origins)
if(unique(st_geometry_type(origins))=="POLYGON"){
if(method=="centroid"){
centers_geom <- st_point_on_surface(origins$geometry)
centers <- origins
centers$geometry <- centers_geom
}else if(method=="pointsalong"){
centers <- surrounding_points(origins,point_dist)
}
}else if(unique(st_geometry_type(origins))=="POINT"){
centers <- origins
}else if(unique(st_geometry_type(origins))=="LINESTRING"){
if(method=="centroid"){
if(verbose){
print("getting the centers of the lines ...")
centers <- lines_center(origins)
}else{
invisible(capture.output(centers <- lines_center(origins)))
}
}else if(method=="ends"){
centers <- lines_extremities(origins)
}else if(method=="pointsalong"){
centers <- lines_points_along(origins,point_dist)
}else{
stop("if origins class is feature collection of linestrings, method must be on of centroid, ends or pointsalong")
}
}else{
stop("origins must be a feature collection of points, polygons or linestrings")
}
if(verbose){
print("snapping the points to the lines (only once)")
}
#snapped_points <- maptools::snapPointsToLines(centers,lines,maxDist = snap_dist, idField="tmpid")
snapped_points <- snapPointsToLines2(centers,lines, idField="tmpid")
#snapped_points <- cbind(snapped_points, centers)
## step 6 building grid
mygrid <- build_grid(grid_shape,list(origins,lines))
if(verbose){
print("starting the network part")
}
ids <- 1:nrow(mygrid)
list_elements <- prepare_elements_netlistw(ids,mygrid,snapped_points,lines,maxdistance)
## step7 iterating over the grid
listvalues <- lapply(1:nrow(mygrid),function(i){
if(verbose){
print(paste("working on quadra : ",i,"/",max(ids),sep=""))
}
elements <- list_elements[[i]]
if(length(elements)==0){
return()
}else {
all_pts <- elements[[1]]
selected_lines <- elements[[2]]
#calculating the elements
values <- network_listw_worker(points = all_pts,
lines = selected_lines,
maxdistance = maxdistance,
direction=direction,
mindist=mindist,
dist_func = vdist_func,
matrice_type = matrice_type,
verbose = verbose,
digits = digits,
tol=tol)
return(values)
}
})
if(verbose){
print("building the final listw object")
}
## step8 combining the results
okvalues <- listvalues[lengths(listvalues) != 0]
nblist <- do.call("c", lapply(okvalues, function(value) {
return(value[[1]])
}))
ordered_nblist <- lapply(as.character(origins$fid), function(x) {
return(nblist[[x]])
})
nbweights <- do.call("c", lapply(okvalues, function(value) {
return(value[[2]])
}))
names(nbweights) <- names(nblist)
ordered_weights <- lapply(as.character(origins$fid), function(x) {
return(nbweights[[x]])
})
# setting the nb attributes
attr(ordered_nblist, "class") <- "nb"
attr(ordered_nblist, "region.id") <- as.character(origins$fid)
if(verbose){
print("finally generating the listw object ...")
}
if (matrice_type != "I"){
# setting the final listw attributes
listw <- spdep::nb2listw(ordered_nblist, glist = ordered_weights,
zero.policy = TRUE, style = matrice_type)
return(listw)
}else{
return(list("nb_list" = ordered_nblist,
"weights" = ordered_weights))
}
}
#' @title Network distance listw (multicore)
#'
#' @description Generate listw object (spdep like) based on network distances with multicore support.
#'
#' @param origins A feature collection of linestrings, points or polygons for
#' which the spatial neighbouring list will be built.
#' @param lines A feature collection of linestrings representing the network
#' @param maxdistance The maximum distance between two observations to consider
#' them as neighbours.
#' @param method A string indicating how the starting points will be built. If
#' 'centroid' is used, then the centre of lines or polygons is used. If
#' 'pointsalong' is used, then points will be placed along polygons' borders
#' or along lines as starting and end points. If 'ends' is used (only for
#' lines) the first and last vertices of lines are used as starting and ending
#' points.
#' @param point_dist A float, defining the distance between points when the
#' method pointsalong is selected.
#' @param snap_dist the maximum distance to snap the start and end points on the
#' network.
#' @param line_weight The weights to use for lines. Default is "length" (the
#' geographical length), but can be the name of a column. The value is
#' considered proportional with the geographical length of the lines.
#' @param mindist The minimum distance between two different observations. It is
#' important for it to be different from 0 when a W style is used.
#' @param direction Indicates a field giving information about authorized
#' travelling direction on lines. if NULL, then all lines can be used in both
#' directions. Must be the name of a column otherwise. The values of the
#' column must be "FT" (From - To), "TF" (To - From) or "Both".
#' @param dist_func Indicates the function to use to convert the distance
#' between observation in spatial weights. Can be 'identity', 'inverse',
#' 'squared inverse' or a function with one parameter x that will be
#' vectorized internally
#' @param matrice_type The type of the weighting scheme. Can be 'B' for Binary,
#' 'W' for row weighted, or 'I' (identity) see the documentation of
#' spdep::nb2listw for details
#' @param grid_shape A vector of length 2 indicating the shape of the grid to
#' use for splitting the dataset. Default is c(1,1), so all the calculation is
#' done in one go. It might be necessary to split it if the dataset is large.
#' @param verbose A Boolean indicating if the function should print its progress
#' @param digits The number of digits to retain in the spatial coordinates (
#' simplification used to reduce risk of topological error)
#' @param tol A float indicating the spatial tolerance when points are added as
#' vertices to lines.
#' @return A listw object (spdep like) if matrice_type is "B" or "W". If
#' matrice_type is I, then a list with a nblist object and a list of weights
#' is returned.
#' @importFrom sf st_length st_point_on_surface
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' future::plan(future::multisession(workers=1))
#' listw <- network_listw.mc(mtl_network,mtl_network,maxdistance=500,
#' method = "centroid", line_weight = "length",
#' dist_func = 'squared inverse', matrice_type='B', grid_shape = c(2,2))
#' ## make sure any open connections are closed afterward
#' if (!inherits(future::plan(), "sequential")) future::plan(future::sequential)
#'}
network_listw.mc <- function(origins,lines,maxdistance, method="centroid", point_dist=NULL, snap_dist=Inf, line_weight = "length", mindist=10, direction=NULL, dist_func = "inverse", matrice_type = "B", grid_shape=c(1,1), verbose = FALSE, digits = 3, tol=0.1){
##adjusting the weights of the lines
lines$line_length <- as.numeric(st_length(lines))
if(line_weight=="length"){
lines$line_weight <- as.numeric(st_length(lines))
}else {
lines$line_weight <- lines[[line_weight]]
}
if(min(lines$line_weight)<=0){
stop("the weights of the lines must be superior to 0")
}
# adjusting the directions of the lines
if(is.null(direction) == FALSE){
lines <- lines_direction(lines,direction)
}
## checking the matrix type
if (matrice_type %in% c("B", "W", "I") == FALSE) {
stop("Matrice type must be B (binary), W (row standardized) or I (identity)")
}
## creating the vectorized distance function
vdist_func <- select_dist_function(dist_func)
## setting the OID : tmpid
lines$tmpid <- 1:nrow(lines)
##generating the starting points
## step5 generating the starting points
if(verbose){
print("generating the starting points")
}
origins$fid <- 1:nrow(origins)
if(unique(st_geometry_type(origins))=="POLYGON"){
if(method=="centroid"){
centers_geom <- st_point_on_surface(origins$geometry)
centers <- origins
centers$geometry <- centers_geom
}else if(method=="pointsalong"){
centers <- surrounding_points(origins,point_dist)
}
}else if(unique(st_geometry_type(origins))=="POINT"){
centers <- origins
}else if(unique(st_geometry_type(origins))=="LINESTRING"){
if(method=="centroid"){
if(verbose){
print("getting the centers of the lines ...")
centers <- lines_center(origins)
}else{
invisible(capture.output(centers <- lines_center(origins)))
}
}else if(method=="ends"){
centers <- lines_extremities(origins)
}else if(method=="pointsalong"){
centers <- lines_points_along(origins,point_dist)
}else{
stop("if origins class is feature collection of linestrings, method must be on of centroid, ends or pointsalong")
}
}else{
stop("origins must be a feature collection of points, polygons or linestrings")
}
if(verbose){
print("snapping the points to the lines (only once)")
}
#snapped_points <- maptools::snapPointsToLines(centers,lines,maxDist = snap_dist, idField="tmpid")
snapped_points <- snapPointsToLines2(centers,lines, idField="tmpid")
#snapped_points <- cbind(snapped_points, centers)
##building grid
grid <- build_grid(grid_shape,list(origins,lines))
if(verbose){
print("starting the network part")
}
all_is <- 1:nrow(grid)
iseq <- list()
cnt <- 0
for(i in 1:grid_shape[[1]]){
start <- cnt*grid_shape[[2]]+1
iseq[[length(iseq)+1]] <- list(cnt+1,all_is[start:(start+grid_shape[[2]]-1)])
cnt<-cnt+1
}
listelements <- future.apply::future_lapply(iseq,function(ii){
elements <- prepare_elements_netlistw(ii[[2]],grid,snapped_points,lines,maxdistance)
return(elements)
}, future.packages = c("sf", "spNetwork"))
listelements <- unlist(listelements,recursive = FALSE)
##iterating on the grid
listvalues <- future.apply::future_lapply(listelements,function(elements){
##step1 : preparing elements
if(is.null(elements)){
return()
}else {
all_pts <- elements[[1]]
selected_lines <- elements[[2]]
#calculating the elements
values <- network_listw_worker(all_pts, selected_lines,
maxdistance, direction=direction, mindist=mindist,
dist_func = vdist_func, matrice_type = matrice_type,
verbose = verbose, digits = digits, tol=tol)
return(values)
}
}, future.packages = c("sf", "spNetwork"))
if(verbose){
print("building the final listw object")
}
okvalues <- listvalues[lengths(listvalues) != 0]
nblist <- do.call("c", lapply(okvalues, function(value) {
return(value[[1]])
}))
ordered_nblist <- lapply(as.character(origins$fid), function(x) {
return(nblist[[x]])
})
nbweights <- do.call("c", lapply(okvalues, function(value) {
return(value[[2]])
}))
names(nbweights) <- names(nblist)
ordered_weights <- lapply(as.character(origins$fid), function(x) {
return(nbweights[[x]])
})
## setting the nb attributes
attr(ordered_nblist, "class") <- "nb"
attr(ordered_nblist, "region.id") <- as.character(origins$fid)
if(verbose){
print("finally generating the listw object ...")
}
if (matrice_type != "I"){
# setting the final listw attributes
listw <- spdep::nb2listw(ordered_nblist, glist = ordered_weights,
zero.policy = TRUE, style = matrice_type)
return(listw)
}else{
return(list("nb_list" = ordered_nblist,
"weights" = ordered_weights))
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.