Merge pull request #6 from RichardPatterson/development

v.0.0.0.9002

.Rbuildignore

@@ -6,3 +6,5 @@
 ^working/examples\.R$
 ^README\.Rmd$
 ^data-raw$
+^doc$
+^Meta$

.gitignore

@@ -1,2 +1,4 @@
 .Rproj.user
 inst/doc
+/doc/
+/Meta/

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: midlines
 Title: Estimate Polygon Midlines
-Version: 0.0.0.9001
+Version: 0.0.0.9002
 Authors@R: 
     person(given = "Richard",
            family = "Patterson",
@@ -12,17 +12,19 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.1
+RoxygenNote: 7.2.3
 Imports: 
     sf,
-    dplyr,
     zoo,
     lwgeom,
-    stats
+    stats,
+    igraph
 Depends: 
     R (>= 2.10)
 Suggests: 
     knitr,
     rmarkdown,
     units
 VignetteBuilder: knitr
+Date: 2024-03-12
+

NEWS.md

@@ -0,0 +1,6 @@
+# midlines 0.0.0.9002
+
+- removed dplyr dependency
+- added igraph dependency
+- replaced wonky bespoke code that grouped intersecting lines with igraph function. This effects midlines_group(), midlines_check(), midlines_dedensify(), midlines_debit(), midlines_smooth(). The new code is ~2.5x fast than the previous version.
+- Vignette added

R/midlines_functions.R

@@ -168,7 +168,7 @@ midlines_clean = function(x, n_removed = 1, border_line = NULL){
 
   }#for loop
 
-  removed_mid_points = dplyr::bind_rows(removed_mid_points)
+  removed_mid_points = do.call("rbind", removed_mid_points)
 
   removed_mid_points$removed_flag = factor(1)
   trimmed_mid_points$removed_flag = factor(0)
@@ -179,7 +179,7 @@ midlines_clean = function(x, n_removed = 1, border_line = NULL){
                        sf::st_drop_geometry(trimmed_points),
                        FUN = unique)
 
-  return = dplyr::inner_join(x, rem_flag, by = "line_id")
+  return = merge(x, rem_flag, by = "line_id")
 
   # Remove line_id var if it wasn't present in input (x)
   if(!line_id_present){
@@ -190,7 +190,7 @@ midlines_clean = function(x, n_removed = 1, border_line = NULL){
 }
 
 
-group_id = line_id = geometry = NULL
+removed_flag = group_id = line_id = geometry = NULL
 
 
 

R/midlines_polish.R

@@ -42,14 +42,16 @@ midlines_debit = function(x, length) {
 midlines_smooth = function(x, width = 3){
 
   dat = midlines_group(x)
-  dat = dplyr::select(dat, geometry)  # stop warning about repeating attributes
-  dat = sf::st_cast(dat,"LINESTRING")
+  dat = subset(dat, select =  geometry)  # stop warning about repeating attributes
+  #dat = sf::st_cast(dat,"LINESTRING")
+  dat = sf::st_sf(geometry = sf::st_cast(sf::st_line_merge(sf::st_union(x)), "LINESTRING"))
+
 
   s = function(x){
     l = length(dat$geometry[[x]])
 
-    a = zoo::rollapply(dat$geometry[[x]][1:(l/2)], width = width, mean )
-    b = zoo::rollapply(dat$geometry[[x]][(l/2+1):l], width = width, mean )
+    a = zoo::rollmean(dat$geometry[[x]][1:(l/2)], k = width)
+    b = zoo::rollmean(dat$geometry[[x]][(l/2+1):l], k = width)
 
     sx = dat$geometry[[x]][1]
     ex = dat$geometry[[x]][(l/2)]
@@ -60,25 +62,20 @@ midlines_smooth = function(x, width = 3){
 
   }
 
-  nrow = nrow(dat)
-
-  smoothed = sf::st_as_sf(sf::st_as_sfc(lapply(1:nrow, s)))
-
-  smoothed = sf::st_as_sf(
-    sf::st_collection_extract(
-      lwgeom::st_split(smoothed$x,
-               sf::st_union(
-                 sf::st_cast(smoothed, "MULTIPOINT"))), type = "LINESTRING"))
-
-  colnames(smoothed)[colnames(smoothed) == "x"] = "geometry"
-  sf::st_geometry(smoothed) <- "geometry"
+  smoothed = sf::st_sfc(lapply(seq_len(nrow(dat)), s))
 
-  smoothed$line_id = 1:nrow(smoothed)
+  smoothed = sf::st_sf(geometry =
+                          sf::st_collection_extract(
+                            lwgeom::st_split(smoothed,
+                                             sf::st_union(
+                                               sf::st_cast(smoothed, "MULTIPOINT"))), type = "LINESTRING"),
+                        crs = sf::st_crs(x))
 
-  sf::st_crs(smoothed) = sf::st_crs(dat)
+  smoothed$line_id = seq_len(nrow(smoothed))
 
   return(smoothed)
-}
+
+  }
 
 
 
@@ -106,8 +103,9 @@ midlines_smooth = function(x, width = 3){
 midlines_dedensify = function(x, density){
 
   x = midlines_group(x)
-  x = dplyr::select(x, geometry)
-  ls = sf::st_cast(x,"LINESTRING")
+  x = subset(x, select = geometry)
+  ls = sf::st_sf(geometry = sf::st_cast(sf::st_line_merge(sf::st_union(x)), "LINESTRING"))
+
   ls$line_id = 1:nrow(ls)
 
   de_densified = sf::st_as_sf(sf::st_line_sample(ls, density = density))

R/midlines_refine.R

@@ -1,4 +1,5 @@
 # An internal function
+# Largely replicates the answer here: https://stackoverflow.com/questions/69175360/is-there-a-way-to-merge-contiguous-multilinestrings-sf-object-in-r
 
 #' Groups line segments into contiguous groups, returns these as multilinestrings.
 #'
@@ -7,48 +8,20 @@
 #' @export
 midlines_group = function(x) {
 
-  lines = sf::st_as_sf(sf::st_cast(sf::st_line_merge(sf::st_union(x)), "LINESTRING"))
+  touches = sf::st_touches(x)
+  graph = igraph::graph_from_adj_list(touches)
 
-  colnames(lines)[colnames(lines) == colnames(lines)] = "geometry" #this only works cos there is one column
-  sf::st_geometry(lines) <- "geometry"
+  groups = igraph::components(graph)$membership
 
-  lines$group_id = NA
-  lines$n_ = 1:nrow(lines)
+  grouped = stats::aggregate(x, by = list(group_id = groups), FUN = unique)
 
-  inter = sf::st_intersects(lines$geometry, lines$geometry)
-  group_index = 1
+  grouped$n_lines = lengths(grouped$line_id)
+  grouped$length = sf::st_length(grouped)
 
-  for(l in 1:nrow(lines)) {
+  #grouped = subset(grouped, select = -c(removed_flag))
+  grouped = grouped[, names(grouped) != "removed_flag"]
 
-    lines$group_id[lines$n_ %in%
-                     inter[[lines$n_[l]]]] = group_index
-
-    lines = lines[order(lines$group_id, na.last = TRUE),]
-
-    if(anyNA(lines$group_id[l+1])) group_index = group_index +1
-    #print(group_index)
-    l = l+1
-  }
-
-  multilines = sf::st_cast(
-    dplyr::summarise(
-      dplyr::group_by(lines, group_id),
-      do_union = FALSE)
-    ,"MULTILINESTRING")
-
-  multilines$n_lines =lengths(lapply(multilines$geometry, unlist))/4
-  multilines$length = sf::st_length(multilines)
-
-  #return(list(multilines,lines))
-  #}
-  grouped_multilinestring = multilines
-  ##############
-
-
-  # this reokaces the variable from the group_lines_a as applying that to that grouped line does something odd.
-  grouped_multilinestring$n_lines = (lengths(lapply(grouped_multilinestring$geometry, unlist)) - 2)/ 2
-
-  return(grouped_multilinestring)
+  return(grouped)
 }
 
 
@@ -107,10 +80,8 @@ midlines_check = function(x, n_removed = NULL, length = NULL, border_line = NULL
   removed = x[x$removed_flag==1,]
   cleaned = x[x$removed_flag==0,]
 
-  #moving this within this cleaning function
   x_multilines = midlines_group(removed)
 
-
   # using n_lines as the number of cycles of removing
   if(!(is.null(n_removed))){
     add_back_groups1 = x_multilines$group_id[x_multilines$n_lines >= n_removed]
@@ -127,7 +98,7 @@ midlines_check = function(x, n_removed = NULL, length = NULL, border_line = NULL
 
   #the first line finds lines touching border and then those removed groups that hit these
   if(!(is.null(border_line))) {
-    #touch_border = cleaned[sf::st_is_within_distance(cleaned, border_line, dist = border_distance, sparse = FALSE),]
+
     touch_border = removed[sf::st_intersects(removed, border_line, sparse = FALSE),]
     add_back_groups3 = x_multilines$group_id[sf::st_intersects(x_multilines, sf::st_union(touch_border), sparse = FALSE)]
   } else {
@@ -136,25 +107,14 @@ midlines_check = function(x, n_removed = NULL, length = NULL, border_line = NULL
 
   add_back_groups = unique(c(add_back_groups1, add_back_groups2, add_back_groups3))
 
-  add_back = x_multilines[x_multilines$group_id %in% add_back_groups,"geometry"]
-
-  # identify the lines to add back from the multilines
-  add_back_index = lengths(sf::st_covered_by(removed, add_back)) != 0
+  add_back_line_ids = unlist(x_multilines$line_id[x_multilines$group_id %in% add_back_groups])
 
-  add_back_lines = removed[add_back_index,]
-  still_removed = removed[!add_back_index,]
+  x$removed_flag2 = x$removed_flag
+  x$removed_flag2[x$line_id %in% add_back_line_ids] = 0
 
-  cleaned$added_flag = factor(0)
-  still_removed$added_flag = factor(0)
-  add_back_lines$added_flag = factor(1)
-
-  cleaned$removed_flag2 = factor(0)
-  add_back_lines$removed_flag2 = factor(0)
-  still_removed$removed_flag2 = factor(1)
-
-  rbind(cleaned, add_back_lines, still_removed)
+  x$added_flag = factor(as.integer(x$removed_flag != x$removed_flag2)) # should probably remove this as its unnecessary
 
+  return(x)
 }
 
 
-

README.Rmd

@@ -50,12 +50,14 @@ library(midlines)
 library(sf)
 library(units)
 
-plot(thames)
+plot(thames$geometry)
+
 ```
 
 To estimate the midline of this stretch of the Thames, use the `midline_draw` function.
 
 ```{r draw}
+
 # create a sf collection for the midline
 m1 = midlines_draw(thames)
 
@@ -115,7 +117,7 @@ bbox_line = st_cast(st_as_sfc(st_bbox(c(xmin = 535070, ymin = 177800, xmax = 542
 
 m1 = midlines_draw(thames, border_line = bbox_line)
 
-plot(thames)
+plot(thames$geometry)
 plot(m1$geometry, col = "BLUE", add = TRUE)
 plot(bbox_line, add = TRUE)
 ```
@@ -142,7 +144,7 @@ bbox_poly = st_as_sfc(st_bbox(c(xmin = 536070, ymin = 180700, xmax = 537800, yma
 # a slightly smaller area to use as a border line
 bbox_line_s = st_cast(st_as_sfc(st_bbox(c(xmin = 536120, ymin = 180760, xmax = 537750, ymax = 181800), crs = st_crs(thames))),"LINESTRING")
 
-plot(thames)
+plot(thames$geometry)
 plot(bbox_poly, add = TRUE)
 ```
 
@@ -151,7 +153,7 @@ plot(bbox_poly, add = TRUE)
 # crop the larger thames polygon to the area we want
 side1 = st_intersection(thames, bbox_poly)
 
-plot(side1)
+plot(side1$geometry)
 plot(bbox_line_s, add = TRUE)
 
 # estimating midline of the original polygon
@@ -163,7 +165,7 @@ The lack of lines is due to the narrowness of the channel, or more precisely, th
 
 ```{r gaps3}
 # some trial and error revealed 8 meters to be the optimal max length between points
-plot(side1)
+plot(side1$geometry)
 plot(bbox_line_s, add = TRUE)
 
 m_s1 = midlines_draw(side1, dfMaxLength = set_units(8,"m"), border_line = bbox_line_s)
@@ -174,7 +176,7 @@ m_s2 = midlines_clean(m_s1, n_removed = 20, border_line = bbox_line_s)
 plot(m_s2$geometry, col = c("BLUE", "RED")[m_s2$removed_flag], add = TRUE)
 ```
 
-Conversely, if the points on the perimeter are too dense and result in a very complicated tessellation or slow computations times then sf::st_line_sample might also be useful.  
+Conversely, if the points on the perimeter are too dense and result in a very complicated tessellation or slow computation times, then sf::st_line_sample might also be useful.  
 
 ## Zig-zagging
 
@@ -227,12 +229,12 @@ bbox_poly = st_as_sfc(st_bbox(c(xmin = 538700, ymin = 180750, xmax = 539800, yma
 # a slightly smaller area to use as a border line
 bbox_line_s = st_cast(st_as_sfc(st_bbox(c(xmin = 538750, ymin = 180800, xmax = 539750, ymax = 181800), crs = st_crs(thames))),"LINESTRING")
 
-plot(thames)
+plot(thames$geometry)
 plot(bbox_poly, add = TRUE)
 
 # crop the larger thames polygon to the area we want
 side2 = st_intersection(thames,bbox_poly)
-plot(side2)
+plot(side2$geometry)
 
 
 plot(bbox_line_s, add = TRUE)
@@ -245,7 +247,7 @@ plot(ml1$geometry, add = TRUE, col = "RED")
 One option would be to use the dfMaxLength option with midlines_draw like we did on the left side channel. 
 
 ```{r smooth}
-plot(side2)
+plot(side2$geometry)
 m_s1a = midlines_draw(side2, border_line = bbox_line_s, dfMaxLength = set_units(8,"m"))
 
 # clean this to remove extraneous lines
@@ -260,7 +262,7 @@ An alternative is to smooth the midline using a rolling average - basically a wr
 # smooth with a rolling average 
 m_s2b = midlines_smooth(m_s1)
 
-plot(side2)
+plot(side2$geometry)
 plot(m_s2b$geometry, add = TRUE, col = "BLUE")
 ```
 
@@ -283,7 +285,7 @@ plot(m2$geometry)
 
 ### De-densifying midlines
 
-In addition to manipulating the density of points on the perimeter of the initial polygon, it is sometimes desirable to manipulate the density of points the form the midline. This results in fewer longer line segments, rather than many shorter line segments. This can be done using `midlines_dedensify` either before or after cleaning. Like with the `midlines_smooth` function, the points where lines join are not affected but between line joins the number of points can be reduced. This function calls sf::st_line_sample so see this for the relevant options.
+In addition to manipulating the density of points on the perimeter of the initial polygon, it is sometimes desirable to manipulate the density of points that form the midline. This results in fewer longer line segments, rather than many shorter line segments. This can be done using `midlines_dedensify` either before or after cleaning. Like with the `midlines_smooth` function, the points where lines join are not affected but between line joins the number of points can be reduced. This function calls sf::st_line_sample so see this for the relevant options.
 
 ## Voronoi polygons
 
@@ -296,6 +298,16 @@ plot(pols)
 
 ```
 
+## Another use case
+
+An alternative use case is to merge two versions of a road network that are measured differently and therefore do not exactly match. Use buffers to creating a polygon around both networks and then estimate the midline of this polygon to give a combined network against which each can be compared. See this [blog](https://richardpatterson.github.io/2023/08/06/network_merge.html) post for more details.
+
+## Citation
+
+If you use this package in your work it would be helpful if you would cite it. A suggested citation is:
+
+_Patterson, R. midlines: Estimate Polygon Midlines. R package version 0.0.0.9002. DOI: 10.17863/CAM.100113_
+
 ## Alternatives
 
 * The very comprehensive R package [`cmgo`](https://github.com/AntoniusGolly/cmgo) does many of the same things as this package and much more. If your aim to to estimate the midline of a river, you should probably start with [`cmgo`](https://github.com/AntoniusGolly/cmgo)