# R/geometrical_functions_sf.R In spNetwork: Spatial Analysis on Network

```# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### helper functions ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#defining some global variables to prevent check error (weird flex but ok)
utils::globalVariables(c("Xs", "Ys", "L1", "lineID"))

#' @title Coordinates to unique character vector
#'
#' @description Generate a character vector based on a coordinates matrix and
#' the maximum number of digits to keep.
#'
#' @param coords A n * 2 matrix representing the coordinates
#' @param digits The number of digits to keep from the coordinates
#' @return A vector character vector of length n
#' @importFrom data.table data.table tstrsplit :=
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
sp_char_index <- function(coords, digits) {
tempdf <- data.frame(Xs = as.character(coords[, 1]),
Ys = as.character(coords[, 2]))
tempdf\$Xs <- ifelse(grepl(".",tempdf\$Xs,fixed = TRUE), tempdf\$Xs, paste(tempdf\$Xs,".0",sep=""))
tempdf\$Ys <- ifelse(grepl(".",tempdf\$Ys,fixed = TRUE), tempdf\$Ys, paste(tempdf\$Ys,".0",sep=""))
tempdt <- data.table(tempdf)

tempdt[,  c("xint", "xdec") := tstrsplit(Xs, ".", fixed = TRUE)]
tempdt[,  c("yint", "ydec") := tstrsplit(Ys, ".", fixed = TRUE)]

X <- paste(tempdt\$xint, substr(tempdt\$xdec, start = 1, stop = digits),
sep = ".")
X <- gsub("NA","0",X,fixed = TRUE)
Y <- paste(tempdt\$yint, substr(tempdt\$ydec, start = 1, stop = digits),
sep = ".")
Y <- gsub("NA","0",Y,fixed = TRUE)
return(paste(X, Y, sep = "_"))
}

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### lines manipulations ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' @title Get lines extremities
#'
#' @description Generate a feature collection of points with the first and last vertex
#'   of each line in a feature collection of linestrings.
#'
#' @param lines A feature collection of linestrings (simple Linestrings)
#' @return A feature collection of points
#' @importFrom sf st_coordinates st_drop_geometry st_crs st_as_sf
#' @importFrom data.table as.data.table .SD .N
#' @export
#' @examples
#' data(mtl_network)
#' points <- lines_extremities(mtl_network)
lines_extremities <- function(lines) {
coords <- st_coordinates(lines)
data <- st_drop_geometry(lines)
remove_it <- FALSE

if(ncol(data) == 1){
data\$jgtmpoid <- 1:nrow(data)
remove_it <- TRUE
}
df <- as.data.table(coords)
pts <- as.data.frame(df[, .SD[c(1,.N)], by=L1])
pts\$pttype <- rep(c("start","end"), nrow(pts)/2)

data <- cbind(data[pts\$L1,], pts[c("pttype","X","Y")])

sf_pts <- st_as_sf(
x = data,
coords = c("X","Y"),
crs = st_crs(lines)
)
sf_pts\$X <- data\$X
sf_pts\$Y <- data\$Y
if(remove_it){
sf_pts\$jgtmpoid <- NULL
}
return(sf_pts)
}

#' @title Remove loops
#'
#' @description Remove from a sf object with linestring type geometries the lines that have the
#' same starting and ending point.
#'
#' @param lines A sf object with linestring type geometries
#' @param digits An integer indicating the number of digits to keep for the
#' spatial coordinates
#' @return A sf object with linestring type geometries
#' @importFrom sf st_coordinates
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
remove_loop_lines <- function(lines, digits){
lines_ext <- lines_extremities(lines)
starts <- subset(lines_ext,lines_ext\$pttype=="start")
ends <- subset(lines_ext,lines_ext\$pttype=="end")
starts\$spoid <- sp_char_index(st_coordinates(starts),digits)
ends\$spoid <- sp_char_index(st_coordinates(ends),digits)
test <- starts\$spoid != ends\$spoid
return(subset(lines,test))
}

#' @title Unify lines direction
#'
#' @description A function to deal with the directions of lines. It ensures
#' that only From-To situation are present by reverting To-From lines. For
#' the lines labelled as To-From, the order of their vertices is reverted.
#'
#' @param lines A sf object with linestring type geometries
#' @param field 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".
#' @return A sf object with linestring type geometries
#' @importFrom sf st_reverse
#' @export
#' @examples
#' data(mtl_network)
#' mtl_network\$length <- as.numeric(sf::st_length(mtl_network))
#' mtl_network\$direction <- "Both"
#' mtl_network[6, "direction"] <- "TF"
#' mtl_network_directed <- lines_direction(mtl_network, "direction")
lines_direction <- function(lines,field){

lines\$tmpjgoid <- 1:nrow(lines)
test <- lines[[field]] == "TF"
rev_sub <- lines[test,]
sub <- lines[!test, ]

reversed <- st_reverse(rev_sub)
combined <- rbind(sub,reversed)
combined <- combined[order(combined\$tmpjgoid),]

return(combined)
}

#' @title Reverse lines
#'
#' @description A function to reverse the order of the vertices of lines
#'
#' @param lines A sf object with linestring type geometries
#' @return A sf object with linestring type geometries
#' @importFrom sf st_coordinates
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
reverse_lines <- function(lines){
return(st_reverse(lines))
}

#' @title Lines coordinates as list
#'
#' @description A function to get the coordinates of some lines as a list of matrices
#'
#' @param lines A sf object with linestring type geometries
#' @return A list of matrices
#' @importFrom sf st_coordinates
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
lines_coordinates_as_list <- function(lines){
line_coords <- st_coordinates(lines)
lineid <- line_coords[,3]
line_coords <- line_coords[,1:2]
line_coords <- split(data.frame(line_coords), f = lineid)
line_coords <- lapply(line_coords,as.matrix)
return(line_coords)
}

# list_coordinates_as_lines_OLD <- function(coord_list, crs){
#
#   lids <- do.call(c,lapply(1:length(coord_list), function(e){rep(e, nrow(coord_list[[e]]))}))
#   new_coords <- data.frame(do.call(rbind, coord_list))
#   new_coords\$lineID <- lids
#   pts <- st_as_sf(new_coords, coords = c(1,2), crs = crs)
#
#   final_lines <- pts %>%
#     group_by(lineID) %>%
#     summarise(do_union = FALSE) %>%
#     st_cast("LINESTRING")
#
#   return(final_lines)
# }

#' @title List of coordinates as lines
#'
#' @description A function to convert a list of matrices to as sf object with linestring geometry type
#'
#' @param coord_list A list of matrices
#' @param crs The CRS to use to create the lines
#' @return A sf object with linestring type geometries
#' @importFrom sf st_coordinates st_cast st_as_sf
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
list_coordinates_as_lines <- function(coord_list, crs){

lids <- do.call(c,lapply(1:length(coord_list), function(e){rep(e, nrow(coord_list[[e]]))}))
new_coords <- data.frame(do.call(rbind, coord_list))
new_coords\$lineID <- lids
colnames(new_coords) <- c('X', 'Y', 'lineID')

final_lines <- sfheaders::sf_linestring(new_coords, x = "X", y = "Y", linestring_id = "lineID", keep = TRUE)
final_lines <- st_as_sf(final_lines)
st_crs(final_lines) <- crs

return(final_lines)
}

#' @title Add vertices to a feature collection of linestrings
#'
#' @description Add vertices (feature collection of points) to their nearest lines
#'   (feature collection of linestrings), may fail if the line geometries are self intersecting.
#'
#' @param lines The feature collection of linestrings to modify
#' @param points The feature collection of points to add to as vertex to the lines
#' @param nearest_lines_idx For each point, the index of the nearest line
#' @param mindist The minimum distance between one point and the extremity of
#'   the line to add the point as a vertex.
#' @return A feature collection of linestrings
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
add_vertices_lines <- function(lines, points, nearest_lines_idx, mindist) {

line_coords <- lines_coordinates_as_list(lines)
pt_coords <- st_coordinates(points)

line_coords,
nearest_lines_idx,
mindist)

final_lines <- list_coordinates_as_lines(new_lines_list, st_crs(lines))
final_lines <- cbind(final_lines, st_drop_geometry(lines))
return(final_lines)
}

# the old version is kept for debug
# lixelize_lines <- function(lines, lx_length, mindist = NULL, verbose = FALSE) {
#     if (is.null(mindist)) {
#         mindist <- lx_length/10
#     }
#     if(verbose){
#         pb <- txtProgressBar(min = 0, max = nrow(lines), style = 3)
#     }
#     newlixels <- lapply(1:nrow(lines), function(i) {
#         if(verbose){
#             setTxtProgressBar(pb, i)
#         }
#         line <- lines[i, ]
#         tot_length <- gLength(line)
#         if (tot_length < lx_length+mindist) {
#             coords <- coordinates(line)
#             return(list(coords[[1]][[1]]))
#         } else {
#             # producing the points to snapp on
#             distances <- seq(lx_length, tot_length, lx_length)
#             if ((tot_length - distances[[length(distances)]]) < mindist) {
#                 distances <- distances[1:(length(distances) - 1)]
#             }
#             points <- t(sapply(distances, function(d) {
#                 return(coordinates(gInterpolate(line, d)))
#             }))
#             points <- data.frame(x = points[, 1], y = points[, 2],
#                                  distance = distances,
#                                  type = "cut")
#             # extracting the original coordinates
#             coords <- SpatialPoints(coordinates(line))
#             xy <- coordinates(coords)
#             points2 <- data.frame(x = xy[, 1],
#                                   y = xy[, 2],
#                                   distance = gProject(line,coords),
#                                   type = "base")
#             if(points2\$distance[nrow(points2)]==0){
#                 points2\$distance[nrow(points2)] <- tot_length
#             }
#             # merging both and sorting
#             allpts <- rbind(points, points2)
#             allpts <- allpts[order(allpts\$distance, allpts\$type), ]
#             # I should remove duplicates here
#             allpts <- allpts[!duplicated(allpts[1:3]),]
#             # and now splitting this motherfucker
#             indices <- c(0, which(allpts\$type == "cut"), nrow(allpts))
#             lixels <- lapply(1:(length(indices) - 1), function(j) {
#                 pts <- allpts[indices[[j]]:indices[[j + 1]], ]
#                 return(as.matrix(pts[, 1:2]))
#             })
#             return(lixels)
#         }
#     })
#     oids <- lapply(1:nrow(lines), function(i) {
#         return(rep(i, length(newlixels[[i]])))
#     })
#     oids <- do.call("c", oids)
#
#     new_lines <- do.call(raster::spLines,unlist(newlixels,recursive = FALSE))
#     df <- lines@data[oids, ]
#     if(class(df) != "dataframe"){
#         df <- data.frame("oid" = df)
#     }
#     new_splines <- SpatialLinesDataFrame(new_lines, df, match.ID = FALSE)
#     raster::crs(new_splines) <- raster::crs(lines)
#     return(new_splines)
#
# }

#' @title Cut lines into lixels
#'
#' @description Cut the lines of a feature collection of linestrings into lixels with a specified minimal
#'   distance may fail if the line geometries are self intersecting.
#'
#' @param lines The sf object with linestring geometry type to modify
#' @param lx_length The length of a lixel
#' @param mindist The minimum length of a lixel. After cut, if the length of the
#'   final lixel is shorter than the minimum distance, then it is added to the
#'   previous lixel. if NULL, then mindist = maxdist/10. Note that the segments
#'   that are already shorter than the minimum distance are not modified.
#' @return An sf object with linestring geometry type
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' lixels <- lixelize_lines(mtl_network,150,50)
#' }
lixelize_lines<- function(lines, lx_length, mindist = NULL) {

if (is.null(mindist)) {
mindist <- lx_length/10
}

if(lx_length <= mindist){
stop("The lixel length (lx_length) must be greater than minimum length (mindist)")
}

coords <- lines_coordinates_as_list(lines)
result <- lixelize_lines_cpp(coords, lx_length, mindist)

final_lines <- list_coordinates_as_lines(result[[1]], st_crs(lines))

df <- st_drop_geometry(lines)[result[[2]]+1, ]

if(inherits(df, "data.frame") == FALSE){
df <- as.data.frame(df)
}

final_lines <- cbind(final_lines, df)

return(final_lines)
}

#'@title Cut lines into lixels (multicore)
#'
#' @description Cut the lines of a feature collection of linestrings into lixels with a specified minimal
#'  distance may fail if the line geometries are self intersecting with
#'  multicore support.
#'
#'@param lines A feature collection of linestrings to convert to lixels
#'@param lx_length The length of a lixel
#'@param mindist The minimum length of a lixel. After cut, if the length of the
#'  final lixel is shorter than the minimum distance, then it is added to the
#'  previous lixel. If NULL, then mindist = maxdist/10
#'@param verbose A Boolean indicating if a progress bar must be displayed
#'@param chunk_size The size of a chunk used for multiprocessing. Default is
#'  100.
#'@return A feature collection of linestrings
#'@export
#'@importFrom utils capture.output
#' @examples
#' \donttest{
#' data(mtl_network)
#' future::plan(future::multisession(workers=1))
#' lixels <- lixelize_lines.mc(mtl_network,150,50)
#' ## make sure any open connections are closed afterward
#' if (!inherits(future::plan(), "sequential")){
#' future::plan(future::sequential)
#' }
#'}
lixelize_lines.mc <- function(lines, lx_length, mindist = NULL, verbose = TRUE, chunk_size = 100) {
chunks <- split(1:nrow(lines), rep(1:ceiling(nrow(lines) / chunk_size),
each = chunk_size, length.out = nrow(lines)))
chunks <- lapply(chunks,function(x){return(lines[x,])})
# step2 : starting the function
iseq <- 1:length(chunks)
if (verbose) {
progressr::with_progress({
p <- progressr::progressor(along = iseq)
values <- future.apply::future_lapply(iseq, function(i) {
p(sprintf("i=%g", i))
chunk_lines <- chunks[[i]]
invisible(capture.output(new_lines <- lixelize_lines(chunk_lines,
lx_length, mindist)))
return(new_lines)
})
})
} else {
values <- future.apply::future_lapply(iseq, function(i) {
chunk_lines <- chunks[[i]]
invisible(capture.output(new_lines <- lixelize_lines(chunk_lines,
lx_length, mindist)))
return(new_lines)
})
}

final_lines <- do.call("rbind", values)
return(final_lines)
}

#' @title LineString to simple Line
#'
#' @description Split the polylines of a feature collection of linestrings in simple
#' segments at each vertex. The values of the columns are duplicated for each segment.
#'
#' @param lines The featue collection of linestrings to modify
#' @return An featue collection of linestrings
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom sf st_crs st_drop_geometry
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' new_lines <- simple_lines(mtl_network)
#' }
simple_lines <- function(lines) {
## extracting the coordinates of the lines
allcoords <- lines_coordinates_as_list(lines)

counts <- sapply(1:length(allcoords), function(i) {
return(nrow(allcoords[[i]]) - 1)
})
oids <- sapply(1:nrow(lines), function(i) {
return(rep(i, counts[[i]]))
})
dim(oids) <- NULL # just to be sure !
if(nrow(lines)>1 & is.list(oids)){
oids <- do.call("c", oids)
}

## using the coordinates to create newlines
new_lines <- lapply(1:length(allcoords), function(i) {
coords <- allcoords[[i]]
segment_lines <- lapply(1:(nrow(coords) - 1), function(i) {
mat <- coords[i:(i + 1), ]
return(mat)
})
return(segment_lines)
})
new_lines <- unlist(new_lines, recursive = FALSE)
final_lines <- list_coordinates_as_lines(new_lines, st_crs(lines))

df2 <- data.frame(st_drop_geometry(lines)[oids,])
names(df2) <- names(st_drop_geometry(lines))
final_lines <- cbind(final_lines, df2)

return(final_lines)
}

#' @title Centre points of lines
#'
#' @description Generate a feature collection of points at the centre of the
#'   lines of a feature collection of linestrings. The length of the lines is
#'   used to determine their centres.
#'
#' @param lines A feature collection of linestrings to use
#' @return A feature collection of points
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' centers <- lines_center(mtl_network)
#' }
lines_center <- function(lines) {

lines\$baseorder <- 1:nrow(lines)

all_lengths <- as.numeric(st_length(lines))
no_length <- subset(lines, all_lengths == 0)
with_length <- subset(lines, all_lengths>0)

## preparing the new points for lines with a length > 0
if(nrow(with_length) > 0){
lines_coords <- lines_coordinates_as_list(with_length)
pts_coords1 <- points_at_lines_centers_cpp(lines_coords)
data1 <- data.frame(st_drop_geometry(with_length))
names(data1) <- names(st_drop_geometry(with_length))
data1\$X <- pts_coords1[,1]
data1\$Y <- pts_coords1[,2]
pts1 <- st_as_sf(data1, coords = c("X","Y"), crs = st_crs(lines))
}else{
pts1 <- NULL
}

## preparing the new points for lines with a length == 0
if(nrow(no_length) > 0){
lines_coords <- lines_coordinates_as_list(no_length)
pts_coords2 <- lapply(lines_coords, function(i){
return(i[1,])
})
pts_coords2 <- do.call(rbind, pts_coords2)
data2 <- data.frame(st_drop_geometry(no_length))
names(data2) <-names(st_drop_geometry(no_length))
data2\$X <- pts_coords2[,1]
data2\$Y <- pts_coords2[,2]
pts2 <- st_as_sf(data2, coords = c("X","Y"), crs = st_crs(lines))
}else{
pts2 <- NULL
}

if(is.null(pts1) == FALSE & is.null(pts2) == FALSE){
all_pts <- rbind(pts1,pts2)
}else if (is.null(pts1) & is.null(pts2) == FALSE){
all_pts <- pts2
} else if (is.null(pts1) == FALSE & is.null(pts2)){
all_pts <- pts1
}else{
stop("impossible to create center points for the provided lines... Check your data")
}

all_pts <- all_pts[order(all_pts\$baseorder),]
all_pts\$baseorder <- NULL

return(all_pts)
}

#' @title Add center vertex to lines
#'
#' @description Add to each feature of a feature collection of lines an additional vertex at its
#' center.
#'
#' @param lines The feature collection of linestrings to use
#' @return A feature collection of points
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
coords <- lines_coordinates_as_list(lines)

final_lines <- list_coordinates_as_lines(new_coords, crs = st_crs(lines))
df <- data.frame(st_drop_geometry(lines))
names(df) <- names(st_drop_geometry(lines))
final_lines <- cbind(final_lines, df)

return(final_lines)
}

#' @title Points along lines
#'
#' @description Generate a feature collection of points along the lines of
#' feature collection of Linestrings.
#'
#' @param lines A feature collection of linestrings to use
#' @param dist The distance between the points along the lines
#' @return A feature collection of points
#' @importFrom utils capture.output
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' new_pts <- lines_points_along(mtl_network,50)
#' }
lines_points_along <- function(lines,dist){

coords <- lines_coordinates_as_list(lines)
new_coords <- points_along_lines_cpp(coords, dist)

df <- data.frame(st_drop_geometry(lines))
deleteit <- FALSE
if(ncol(df) == 1){
df\$tmpjgid <- 0
deleteit <- TRUE
}
df <- df[new_coords[,3]+1,]
if(deleteit){
df\$tmpjgid <- NULL
}
df\$X <- new_coords[,1]
df\$Y <- new_coords[,2]
all_pts <- st_as_sf(df, coords = c("X","Y"), crs = st_crs(lines))

return(all_pts)
}

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### functions on polygons ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' @title Points along polygon boundary
#'
#' @description Generate a feature collection of points by placing points along the border of
#' polygons of a feature collection.
#'
#' @param polygons A feature collection of polygons
#' @param dist The distance between the points
#' @return A feature collection of points representing the points arrond the polygond
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
surrounding_points <- function(polygons,dist){
boundaries <- st_boundary(polygons)
pts <- lines_points_along(boundaries,dist)
return(pts)
}

# this is the previous version of the function, kept for debuging
# lines_points_along <- function(lines,dist){
#     lenghts <- gLength(lines, byid = TRUE)
#     list_pts <- lapply(1:nrow(lines),function(i){
#         line <- lines[i,]
#         line_lenght <- lenghts[i]
#         distances <- seq(0,line_lenght,dist)
#         pts <- gInterpolate(line,distances)
#         return(pts)
#     })
#     oids <- lapply(1:length(list_pts),function(i){rep(i,length(list_pts[[i]]))})
#     oids <- do.call("c",oids)
#     all_pts <- do.call(rbind,list_pts)
#     # adding a useless column to avoid a bug when lines has only one column
#     lines\$tmpOID <- 1:nrow(lines)
#     data <- lines@data[oids,]
#     all_pts <- sp::SpatialPointsDataFrame(all_pts,data)
#     all_pts\$tmpOID <- NULL
#     return(all_pts)
# }

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### gridding function ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' @title sf geometry bbox
#'
#' @description Generate polygon as the bounding box of a feature collection
#'
#' @param x A feature collection
#' @return A feature collection of polygons
#' @importFrom sf st_as_sfc st_bbox
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
st_bbox_geom <- function(x){
vec <- st_bbox(x)
poly <- st_as_sf(st_as_sfc(vec),crs = st_crs(x))
poly\$oid <- 1
return(poly)
}

#' @title Spatial grid
#'
#' @description Generate a grid of a specified shape in the bbox of a Spatial object.
#'
#' @param grid_shape A numeric vector of length 2 indicating the number
#' of rows and the numbers of columns of the grid
#' @param spatial A list of spatial feature collections objects (package sf)
#' @return A feature collection of polygons representing the grid
#' @importFrom sf st_make_grid
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
build_grid <- function(grid_shape, spatial) {
boxes <- lapply(spatial,function(i){matrix(st_bbox(i),nrow = 2, byrow = FALSE)})
boxes <- do.call(cbind,boxes)
v1 <- as.numeric(c(min(boxes[1,]),max(boxes[1,])))
v2 <- as.numeric(c(min(boxes[2,]),max(boxes[2,])))
box <- rbind(v1,v2)
vec <- c(box[1,1], box[2,1], box[1,2], box[2,2])
names(vec) <- c("xmin", "ymin", "xmax", "ymax")
class(vec) <- "bbox"
if(prod(grid_shape)==1){
poly <- st_as_sf(st_as_sfc(vec),crs = st_crs(spatial[[1]]))
poly\$oid <- 1
return(poly)
}else{
## step1 : creating the grid
grid <- st_as_sf(st_make_grid(vec, n = grid_shape,
crs = st_crs(spatial[[1]]),
what = "polygons"))
grid\$oid <- 1:nrow(grid)
return(grid)
}
}

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### snapping function ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' @title Nearest point on segment
#'
#' @description Find the nearest projected point on a segment (from maptools)
#'
#' @param s The coordinates of the segment
#' @param p The coordinates of the point
#'
#' @return A numeric vector with the coordinates of the projected point
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
nearestPointOnSegment <- function(s, p){
ap <- c(p[1] - s[1,1], p[2] - s[1,2])
ab <- c(s[2,1] - s[1,1], s[2,2] - s[1,2])
t <- sum(ap*ab) / sum(ab*ab)
t <- ifelse(t<0,0,ifelse(t>1,1,t))
# when start and end of segment are identical t is NA
t <- ifelse(is.na(t), 0, t)
x <- s[1,1] + ab[1] * t
y <- s[1,2] + ab[2] * t
# Return nearest point and distance
result <- c(x, y, sqrt((x-p[1])^2 + (y-p[2])^2))
names(result) <- c("X","Y","distance")
return(result)
}

#' @title Nearest point on Line
#'
#' @description Find the nearest projected point on a LineString (from maptools)
#'
#' @param coordsLine The coordinates of the line (matrix)
#' @param coordsPoint The coordinates of the point
#'
#' @return A numeric vector with the coordinates of the projected point
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
nearestPointOnLine <- function(coordsLine, coordsPoint){
nearest_points <- vapply(2:nrow(coordsLine), function(x){
nearestPointOnSegment(coordsLine[(x-1):x,], coordsPoint)
}, FUN.VALUE=c(0,0,0))

# Return coordinates of the nearest point on this line
return(nearest_points[1:2, which.min(nearest_points[3,])])

}

# NOTE: this code is kept if bugs are encountered in the new version using
# @title Nearest feature
#
# @description Find the nearest feature from set X in set Y
#
# @param x A feature collection of points
# @param y A feature collection of points
#
# @return A numeric vector with the index of the nearest features
# @keywords internal
# @examples
# #This is an internal function, no example provided
# nearest <- function(x,y){
#     sf1 <- sf::st_as_sf(x)
#     sf2 <- sf::st_as_sf(y)
#     sf::st_crs(sf1) <- sf::st_crs(sf2)
#     best <- sf::st_nearest_feature(sf1,sf2)
#     return(best)
# }

#' @title Nearest line for points
#' @description Find for each point its nearest LineString
#' @param points A feature collection of points
#' @param lines A feature collection of linestrings
#' @param snap_dist A distance (float) given to find for each point its
#' nearest line in a spatial index. A too big value will produce
#' unnecessary distance calculations and a too short value will lead to
#' more iterations to find neighbours. In extrem cases, a too short value
#' could lead to points not associated with lines (index = -1).
#' @param max_iter An integer indicating how many iteration the search
#' algorithm must perform in the spatial index to find lines close to a
#' point. At each iteration, the snap_dist is doubled to find candidates.
#' @keywords internal
#' @examples
#' # this is an internal function, no example provided
nearest_lines <- function(points, lines, snap_dist = 300, max_iter = 10){

# getting the coordinates of the lines
list_lines <- lines_coordinates_as_list(lines)

# getting the coordinates of the points
coords <- st_coordinates(points)

# getting the indexes
idx <- find_nearest_object_in_line_rtree(coords, list_lines, snap_dist, max_iter)

# adding 1 to match with c++ indexing
return(idx+1)
}

# snapPointsToLines2_OLD <- function(points, lines ,idField = NA, snap_dist = 300, max_iter = 10) {
#
#     #nearest_line_index <- nearest(points,lines)
#     if(is.na(idField)){
#       lines\$tmpjgid <- 1:nrow(lines)
#       idField <- "tmpjgid"
#     }
#
#     nearest_line_index <- nearest_lines(points, lines, snap_dist, max_iter)
#     coordsLines <- lines_coordinates_as_list(lines)
#     coordsPoints <- st_coordinates(points)
#
#     # Get coordinates of nearest points lying on nearest lines
#     mNewCoords <- vapply(1:nrow(points),function(x){
#             return(nearestPointOnLine(coordsLines[[nearest_line_index[x]]],
#                                coordsPoints[x,]))}, FUN.VALUE=c(0,0))
#
#     # Recover lines' Ids (If no id field has been specified, take the sp-lines id)
#     nearest_line_id <- lines[[idField]][nearest_line_index]
#
#     # Create data frame and sp points
#     df <- data.frame(st_drop_geometry(points))
#
#     df\$nearest_line_id <- nearest_line_id
#     mNewCoords <- t(mNewCoords)
#     df\$Xx <- mNewCoords[,1]
#     df\$Yy <- mNewCoords[,2]
#     final_points <- st_as_sf(df, coords = c("Xx","Yy"), crs = st_crs(lines))
#     final_points\$Xx <- NULL
#     final_points\$Yy <- NULL
#
#     return(final_points)
#
# }

#' @title Snap points to lines
#'
#' @description Snap points to their nearest lines (edited from maptools)
#'
#' @param points A feature collection of points
#' @param lines A feature collection of linestrings
#' @param idField The name of the column to use as index for the lines
#' @param snap_dist A distance (float) given to find for each point its
#' nearest line in a spatial index. A too big value will produce
#' unnecessary distance calculations and a too short value will lead to
#' more iterations to find neighbours. In extrem cases, a too short value
#' could lead to points not associated with lines (index = -1).
#' @param max_iter An integer indicating how many iteration the search
#' algorithm must perform in the spatial index to find lines close to a
#' point. At each iteration, the snap_dist is doubled to find candidates.
#'
#' @return A feature collection of points with the projected geometries
#' @keywords internal
#' @importFrom methods slot
#' @export
#' @examples
#' data(mtl_network)
#' data(bike_accidents)
#' mtl_network\$LineID <- 1:nrow(mtl_network)
#' # snapping point to lines
#' snapped_points <- snapPointsToLines2(bike_accidents,
#'     mtl_network,
#'     "LineID"
#' )
snapPointsToLines2 <- function(points, lines ,idField = NA, snap_dist = 300, max_iter = 100) {

#nearest_line_index <- nearest(points,lines)
if(is.na(idField)){
lines\$tmpjgid <- 1:nrow(lines)
idField <- "tmpjgid"
}

#nearest_line_index <- nearest_lines(points, lines, snap_dist, max_iter)

# we start by finding the nearest line for each point
nearest_line_index <- sf::st_nearest_feature(points, lines)
lines2 <- lines[nearest_line_index,]

# we can then project the points on theses lines
proj_lines <- sf::st_nearest_points(points, lines2, pairwise = TRUE, tolerance = 1e-9)
proj_pts <- st_cast(proj_lines, "POINT")
final_points <- st_as_sf(proj_pts[seq(2,length(proj_pts),2)])

# we can then add the ids of the lines
final_points\$nearest_line_id <- lines2[[idField]]

# and we can create the final object
final_points <- cbind(final_points, st_drop_geometry(points))
st_geometry(final_points) <- 'geometry'

return(final_points)

}

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### Network Cleaning Functions ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#' @title Heal edges
#'
#' @description Merge Lines if they form a longer linestring without external intersections (experimental)
#'
#' @param lines A feature collection of linestrings
#' @param digits An integer indicating the number of digits to keep in coordinates
#' @param verbose A boolean indicating if a progress bar should be displayed
#'
#' @return A feature collection of linestrings with the eventually merged geometries. Note
#' that if lines are merged, only the attributes of the first line are preserved
#' @keywords internal
#' @importFrom sf st_line_merge st_sfc st_cast st_geometry st_geometry<-
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @examples
#' #This is an internal function, no example provided
heal_edges <- function(lines,digits = 3, verbose = TRUE){
## healing step
geoms <- st_geometry(lines)
lines <- st_drop_geometry(lines)
sf::st_geometry(lines) <- geoms

#first getting the coordinates of each point
lines\$tmpOID <- 1:nrow(lines)
coords <- lines_extremities(lines)
coords\$spindex <- sp_char_index(st_coordinates(coords),digits)

#then counting for each point how many times it appears
counts <- table(coords\$spindex)
countdf <- data.frame(
count = as.vector(counts),
spindex = names(counts)
)

#if a point appears exactly twice, it might require a healing
cases <- subset(countdf, countdf\$count==2)
consid_points <- subset(coords, coords\$spindex %in% cases\$spindex)

#adding the start and end sp index to each line
startpts <- data.frame(st_drop_geometry(subset(coords, coords\$pttype=="start")))
endpts <- data.frame(st_drop_geometry(subset(coords, coords\$pttype=="end")))
names(startpts) <- names(st_drop_geometry(coords))
names(endpts) <- names(startpts)

tempDT <- as.data.table(lines)
lines\$startidx <- setDT(tempDT)[startpts, on = "tmpOID", "startidx" := startpts\$spindex][]\$startidx
lines\$endidx <- setDT(tempDT)[endpts, on = "tmpOID", "endidx" := endpts\$spindex][]\$endidx

test <- (lines\$startidx %in% consid_points\$spindex | lines\$endidx %in% consid_points\$spindex)
consid_lines <- subset(lines, test)
keeped_lines <- subset(lines, test == FALSE)

consid_lines\$tmpOID <- as.character(consid_lines\$tmpOID)
# generating a neighbouring list from the extremies

neighbouring <- lapply(1:nrow(consid_lines), function(i){
line <- consid_lines[i,]
val1 <- line\$startidx
val2 <- line\$endidx
if(val1 %in% consid_points\$spindex){
test1 <- (consid_lines\$startidx == val1 | consid_lines\$endidx == val1) &
((consid_lines\$startidx == val1 & consid_lines\$endidx == val2)==FALSE
)
}else{
test1 <- rep(FALSE, nrow(consid_lines))
}
if(val2 %in% consid_points\$spindex){
test2 <- (consid_lines\$startidx == val2 | consid_lines\$endidx == val2) &
((consid_lines\$startidx == val1 & consid_lines\$endidx == val2)==FALSE
)
}else{
test2 <- rep(FALSE, nrow(consid_lines))
}
neighbours <- subset(consid_lines, test1 | test2)

codes <- neighbours\$tmpOID
return(codes)
})
names(neighbouring) <- consid_lines\$tmpOID

## ok now we must find the routes that we should merge
merged <- c()
if(verbose){
pb <- txtProgressBar(min = 0, max = nrow(consid_lines), style = 3)
}
merge_with <- lapply(1:nrow(consid_lines), function(i){
line <- consid_lines[i,]
if(verbose){
setTxtProgressBar(pb, i)
}
all_neighbours <- c()
next_check <- line\$tmpOID
while (length(next_check) > 0){
neighbours <- do.call(c,lapply(next_check, function(x){
neighbouring[[x]]
}))
neighbours <- neighbours[(neighbours %in% merged) == FALSE]
next_check <- neighbours
merged <<- c(merged,neighbours)
all_neighbours <- c(all_neighbours,neighbours)
}
return(all_neighbours)
})
merge_with <- Filter(function(x) length(x) > 0, merge_with)

merged_lines <- st_sfc(do.call(rbind,lapply(merge_with, function(x){
sub <- subset(consid_lines, consid_lines\$tmpOID %in% x)
return(st_line_merge(st_cast(st_geometry(sub),"MULTILINESTRING", ids = rep(1, nrow(sub)))))
})))

df <- do.call(rbind,lapply(merge_with, function(x){
sub <- subset(consid_lines, consid_lines\$tmpOID %in% x)
return(sub[1,])
}))

df\$geometry <- merged_lines
df\$startidx <- NULL
df\$endidx <- NULL
keeped_lines\$startidx <- NULL
keeped_lines\$endidx <- NULL
st_crs(df) <- st_crs(keeped_lines)
new_lines <- rbind(keeped_lines, df)
new_lines\$tmpOID <- NULL
return(new_lines)
}

#' @title Remove mirror edges
#'
#' @description Keep unique edges based on start and end point
#'
#' @param lines A feature collection of linestrings
#' @param keep_shortest A boolean, if TRUE, then the shortest line is keeped if
#' several lines have the same starting point and ending point. if FALSE, then the
#' longest line is keeped.
#' @param digits An integer indicating the number of digits to keep in coordinates
#'
#' @return A feature collection of linestrings with the mirror edges removed
#' @keywords internal
#' @importFrom sf st_length
#' @importFrom data.table as.data.table setDT :=
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @examples
#' #This is an internal function, no example provided
remove_mirror_edges <- function(lines, keep_shortest = TRUE, digits = 3, verbose = TRUE){

# step1 : get start and end coordinates of each line
lines\$tmpOID <- 1:nrow(lines)
coords <- lines_extremities(lines)
coords\$spindex <- sp_char_index(st_coordinates(coords),digits)

starts <- subset(coords, coords\$pttype == "start")
starts\$geometry <- NULL
ends <- subset(coords, coords\$pttype == "end")
ends\$geometry <- NULL
# merging with data.table
dfline <- lines
dfline\$geometry <- NULL

tempDT <- data.table::as.data.table(dfline)
lines\$startidx <- setDT(tempDT)[starts, on = "tmpOID", "startidx" := starts\$spindex][]\$startidx
lines\$endidx <- setDT(tempDT)[ends, on = "tmpOID", "endidx" := ends\$spindex][]\$endidx

# create two sp_index (one start->end and end->start)
lines\$idx1 <- paste(lines\$startidx,lines\$endidx,sep=" - ")
lines\$idx2 <- paste(lines\$endidx,lines\$startidx,sep=" - ")

# step2 : for each line, find if another has the same paire of points (or reversed)
# if it is the case, flag the longest or shortest for a removal after
if(verbose){
pb <- txtProgressBar(min = 0, max = nrow(lines), style = 3)
}
to_remove <- sapply(1:nrow(lines), function(i){
line <- lines[i,]
idx1 <- line\$idx1
if(verbose){
setTxtProgressBar(pb, i)
}
test <- (lines\$idx1 == idx1 | lines\$idx2 == idx1)
if (sum(test)>1){
sub <- subset(lines, test)
lengths <- as.numeric(st_length(sub))
# if all the lengths are equal
if(length(unique(lengths))==1){
if(sum(line\$tmpOID > sub\$tmpOID)==0){
return(TRUE)
}else{
return(FALSE)
}
}
# if some lengths are longer
l <- as.numeric(st_length(line))
longest <- max(lengths)
if(keep_shortest){
if(longest > l){
return(FALSE)
}else{
return(TRUE)
}
}else{
if(longest < l){
return(FALSE)
}else{
return(TRUE)
}
}
}else{
return(FALSE)
}
})
keeped <- subset(lines, to_remove == FALSE)
return(keeped)
}

#' @title Simplify a network
#'
#' @description Simplify a network by applying two corrections: Healing edges and
#' Removing mirror edges (experimental).
#'
#' @details Healing is the operation to merge two connected linestring if the are
#' intersecting at one extremity and do not intersect any other linestring. It helps
#' to reduce the complexity of the network and thus can reduce calculation time.
#' Removing mirror edges is the operation to remove edges that have the same
#' extremities. If two edges start at the same point and end at the same point,
#' they do not add information in the network and one can be removed to simplify
#' the network. One can decide to keep the longest of the two edges or the shortest.
#' NOTE: the edge healing does not consider lines directions currently!
#'
#' @param lines A feature collection of linestrings
#' @param digits An integer indicating the number of digits to keep in coordinates
#' @param heal A boolean indicating if the healing operation must be performed
#' @param mirror A boolean indicating if the mirror edges must be removed
#' @param keep_shortest A boolean, if TRUE, then the shortest line is kept from
#' mirror edges. if FALSE, then the longest line is kept.
#' @param verbose A boolean indicating if messages and a progress bar should be displayed
#' @return A feature collection of linestrings
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' edited_lines <- simplify_network(mtl_network, digits = 3, verbose = FALSE)
#' }
simplify_network <- function(lines, digits = 3, heal = TRUE, mirror = TRUE, keep_shortest = TRUE, verbose = TRUE){ # nocov start

if(heal){
if(verbose){
print("healing the connected edges...")
}
lines <- heal_edges(lines, digits, verbose)
}

if(mirror){
if(verbose){
print("removing mirror edges")
}
lines <- remove_mirror_edges(lines, keep_shortest, digits, verbose)
}

return(lines)
} # nocov end

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### Development ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# # # this is the old function kept if bugs are found
# split_at_vertices <- function(line, points, i, mindist) {
#     # extract coordinates
#     line_coords <- coordinates(line)[[1]][[1]]
#     original_distances <- sapply(1:nrow(line_coords),function(i){
#         if (i==0){
#             return(0)
#         }else{
#             return(sqrt(sum((line_coords[i,]-line_coords[i-1,])**2)))
#         }
#     })
#     line_coords <- cbind(line_coords, cumsum(original_distances),0)
#     tot_lengths <- max(line_coords[,3])
#     # calculate lengths
#     pt_coords <- coordinates(points)
#     lengths <- gProject(line, points)
#     pt_coords <- cbind(pt_coords, lengths)
#     pt_coords <- pt_coords[(lengths>mindist & lengths<(tot_lengths-mindist)),]
#     if(is.null(nrow(pt_coords))){
#         pt_coords <- c(pt_coords,1)
#     }else{
#         pt_coords <- cbind(pt_coords,1)
#     }
#
#     all_coords <- rbind(line_coords,pt_coords)
#     # reorder the coordinate matrix
#     ord_coords <- all_coords[order(all_coords[,3]), ]
#
#     #split on new coordinates
#     ruptidx <- unique(c(1,(1:nrow(ord_coords))[ord_coords[,4] == 1],nrow(ord_coords)))
#     final_coords <- lapply(1:(length(ruptidx)-1), function(j){
#         els <- ord_coords[ruptidx[[j]]:ruptidx[[j+1]],1:2]
#     })
#
#     return(final_coords)
# }

# # this is the previous function, kept for debug
# split_lines_at_vertex2 <- function(lines, points, nearest_lines_idx, mindist) {
#     new_lines_list <- lapply(1:nrow(lines), function(i) {
#         line <- lines[i, ]
#         testpts <- nearest_lines_idx == i
#         if (any(testpts)) {
#             okpts <- subset(points,testpts)
#             newline <- split_at_vertices(line, okpts, i, mindist)
#             return(newline)
#         } else {
#             sline <- list(coordinates(line)[[1]][[1]])
#             return(sline)
#         }
#     })
#     final_lines <- do.call(raster::spLines,unlist(new_lines_list, recursive = FALSE))
#     idxs <- do.call(c,lapply(1:length(new_lines_list), function(j){
#         el <- new_lines_list[[j]]
#         if (class(el) == "list"){
#             return(rep(j,times = length(el)))
#         }else{
#             return(j)
#         }
#     }))
#     final_lines <- SpatialLinesDataFrame(final_lines,
#                                          lines@data[idxs,],match.ID = FALSE)
#     raster::crs(final_lines) <- raster::crs(lines)
#     return(final_lines)
# }

#' @title Split lines at vertices in a feature collection of linestrings
#'
#' @description Split lines (feature collection of linestrings) at their nearest vertices
#' (feature collection of points), may fail if the line geometries are self intersecting.
#'
#' @param lines The feature collection of linestrings to split
#' @param points The feature collection of points to add to as vertex to the lines
#' @param nearest_lines_idx For each point, the index of the nearest line
#' @param mindist The minimum distance between one point and the extremity of
#'   the line to add the point as a vertex.
#' @return A feature collection of linestrings
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
#' @examples
#' \donttest{
#' data(mtl_network)
#' data(bike_accidents)
#' # aggregating points within a 5 metres radius
#' bike_accidents\$weight <- 1
#' agg_points <- aggregate_points(bike_accidents, 5)
#' mtl_network\$LineID <- 1:nrow(mtl_network)
#' # snapping point to lines
#' snapped_points <- snapPointsToLines2(agg_points,
#'     mtl_network,
#'     "LineID"
#' )
#' # splitting lines
#' new_lines <- split_lines_at_vertex(mtl_network, snapped_points,
#'     snapped_points\$nearest_line_id, 1)
#' }
split_lines_at_vertex <- function(lines, points, nearest_lines_idx, mindist) {

coords <- st_coordinates(points)
# step1 : remove points that are not far enough from the extremities
exts <- lines_extremities(lines)

coords2 <- st_coordinates(exts)

#min_dists <- FNN::knnx.dist(data = coords2, query = coords,k = 1)
min_dists <- dbscan::kNN(x = coords2, query = coords,k = 1)\$dist
coords <- subset(coords, min_dists > mindist)
nearest_lines_idx <- nearest_lines_idx[min_dists > mindist]
# no need to split here, great !
if(nrow(coords) == 0){
return(lines)
}else{
#lines_coords <- unlist(sp::coordinates(lines), recursive = FALSE)
lines_coords <- lines_coordinates_as_list(lines)
elements <- split_lines_at_points_cpp(coords, lines_coords, nearest_lines_idx, mindist)
new_lines_list <- elements[[1]]
idxs <- elements[[2]]

#final_lines <- do.call(raster::spLines,new_lines_list)
final_lines <- list_coordinates_as_lines(new_lines_list, crs = st_crs(lines))
final_lines <- cbind(final_lines, st_drop_geometry(lines[idxs,]))
return(final_lines)
}
}

#' @title Cut lines at a specified distance
#'
#' @description Cut lines in a feature collection of linestrings at a specified distance from the
#' begining of the lines.
#'
#' @param lines The feature collection of linestrings to cut
#' @param dists A vector of distances, if only one value is given,
#' each line will be cut at that distance.
#' @return A feature collection of linestrings
#' @keywords internal
#' @examples
#' # This is an interal function, no example provided
cut_lines_at_distance <- function(lines, dists){

# step 1 : create a list of coordinates
coord_lists <- lines_coordinates_as_list(lines)

new_coords <- cut_lines_at_distances_cpp(coord_lists, dists)

final_lines <- list_coordinates_as_lines(new_coords, crs = st_crs(lines))
final_lines <- cbind(final_lines, st_drop_geometry(lines))

return(final_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.