Connect Linestrings

Question

How can i connect two linestrings?

It is possible to lay a very slim buffer around the lines and then connect them like so:

one_line <- lines %>% 
  st_buffer(0.05) %>% 
  st_intersection() %>% 
  st_union() %>% 
  st_cast('LINESTRING')

There are 2 problems with this:

a) below is a very small subset of my data containing one such disconnected line segment - if i use the above method on the small part it forms a complete polygon which, when converted to a linestring just makes a very narrow loop

b) if i use the whole data set it kind of works but creates lines at the approximate distance of the buffer around my original line. See picture below:

Blue & red are the edge lines while black would be the original.

I thought to simply average them out but when i convert the 2 lines to coordinates (st_coordinates()), the resulting tables have different lengths and are not in order.

I looked around but did not really find any useful answers.

Here is a dput of the geometry data:

lines <- structure(list(structure(list(structure(c(2880, 2880.92, 2881.72, 
                                                2882.47, 2883.17, 2883.84, 2884.5, 2894.05, 2894.69, 2895.29393034826, 
                                                340255.362641509, 340257.22, 340259.03, 340260.85, 340262.69, 
                                                340264.55, 340266.4, 340293.7, 340295.61, 340297.500995024), .Dim = c(10L, 
                                                                                                                      2L)), structure(c(2907.22402724177, 2914.21353757771, 340330.886392736, 
                                                                                                                                        340350.2), .Dim = c(2L, 2L))), class = c("XY", "MULTILINESTRING", 
                                                                                                                                                                                 "sfg")), structure(c(2895.3, 2896.82, 2897.26, 2897.72, 2907.2, 
                                                                                                                                                                                                      340297.52, 340302.26, 340303.58, 340304.89, 340330.82), .Dim = c(5L, 
                                                                                                                                                                                                                                                                       2L), class = c("XY", "LINESTRING", "sfg"))), n_empty = 0L, crs = structure(list(
                                                                                                                                                                                                                                                                         input = "EPSG:31256", wkt = "PROJCRS[\"MGI / Austria GK East\",\n    BASEGEOGCRS[\"MGI\",\n        DATUM[\"Militar-Geographische Institut\",\n            ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4312]],\n    CONVERSION[\"Austria Gauss-Kruger East\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",16.3333333333333,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",1,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-5000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Austria east of 14Â°50'E of Greenwich (32Â°30'E of Ferro).\"],\n        BBOX[46.56,14.83,49.02,17.17]],\n    ID[\"EPSG\",31256]]"), class = "crs"), idx = structure(c(1, 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             2, 1, 1), .Dim = c(2L, 2L)), class = c("sfc_GEOMETRY", "sfc"), precision = 0, bbox = structure(c(xmin = 2880, 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              ymin = 340255.362641509, xmax = 2914.21353757771, ymax = 340350.2
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             ), class = "bbox"), classes = c("MULTILINESTRING", "LINESTRING"
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             ))

Which lines do you want to connect, and how do you want them to be connected? — mrhellmann, Mar 08 '22 at 14:05
well it is almost impossible to notice but in the dput sample but if you plot it like so: `x11(); plot(lines[1]); plot(lines[2], col="red", add=T)` you can see that the linesegments do not actually touch (the difference is a miniscule 0.02m) but with `zoom`s `zm()` function you can zoom in near enough to see that there is a gap. i would like to connect those two lines to get a single linestring. - right now i am trying to achive this by finding the closest end-/startpoint (st_endpoint/st_startpoints) — D.J, Mar 08 '22 at 14:10
Hi D.J I understand in the comments below that your dataset is a bit more complex. Would it be possible to modify your post by adding a more representative sample of your dataset in order to take into account all the cases you are facing (e.g. zig-zag lines)? Cheers. — lovalery, Mar 11 '22 at 17:23
@lovalery sadly no since they are work files which i am not at liberty to share (small parts are ok but the whole file is a no go) else i would have shared a full sample of the geometry. are you interested in the problem for fun or do you habe a similar problem? — D.J, Mar 11 '22 at 19:53
@D.J thanks for your feedback. I understand. No problem. It was just for fun and to help you in case you need an alternative solution. Cheers. — lovalery, Mar 11 '22 at 20:18

score 1 · Answer 1 · answered Mar 08 '22 at 14:59

1

For the above example, you could cast to MULTIPOINT, then union, and cast to LINESTRING.

``` r
library(tidyverse) #overkill, but easier
library(sf)
library(patchwork) #to plot side-by-side

# load data from above 
# lines <- 

single_line <- lines %>%
  st_as_sf() %>%
  st_cast('MULTIPOINT') %>%
  st_union() %>%
  st_cast('LINESTRING')

head(single_line)
#> Geometry set for 1 feature 
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 2880 ymin: 340255.4 xmax: 2914.214 ymax: 340350.2
#> Projected CRS: MGI / Austria GK East
#> LINESTRING (2880 340255.4, 2880.92 340257.2, 28...

p1 <- ggplot() +
  geom_sf(data = st_as_sf(lines), col = c('red', 'blue')) + ggtitle('lines')
          
p2 <- ggplot() + geom_sf(data = single_line, col = 'black') + ggtitle('lines cast & unioned')

p3 <- p1 + p2

p3

^{Created on 2022-03-08 by the reprex package (v0.3.0)}

answered Mar 08 '22 at 14:59

mrhellmann

5,069
11
38

1

i will check this when i get back to work, casting the line to multipoint seems like a great idea! thank you! – D.J Mar 08 '22 at 15:09
i tested this and for my small example it works perfectly. for my larger data set the points are not ordered and i get a zig-zag line instead of the correct geometry. – D.J Mar 09 '22 at 05:54
1

@D.J I thought that might happen. I'm pretty sure there are ways around that. Maybe post a new question with some data that causes the problem? I'll take another look & it'll likely get a few other `sf` watchers to do the same. – mrhellmann Mar 09 '22 at 11:01
There have to be. Tbh i did not think that such a simple problem would turn out to be so complex. I think i will post another question with the data of a full line - although duplicate questions are discouraged. – D.J Mar 09 '22 at 11:09
I found a solution - it is not perfect because it still needs some manual labor but for simple cases it works well. – D.J Mar 09 '22 at 14:56
1

@D.J Good to hear. If there's too much manual labor, there are probably workarounds with buffers, line sampling, sfheaders, and/or maybe networks. – mrhellmann Mar 09 '22 at 15:03
there is a "workaround". i spent most of the day writing some iterator that goes through all the line segments and pairs them. since the differences are in the submeter area i can filter accordingly. but it is ugly, convoluted and barely tested so i didnt include it here :) anyway your idea really helped me out so thanks again! (used in in the answer) – D.J Mar 09 '22 at 15:06

score 0 · Answer 2 · answered Mar 09 '22 at 15:02

This is an incomplete answer because it still requires some manual input but it can be generalized by implementing a few rules. The idea is to find the places where the line is broken and then go to the start/endpoints of the closest line. Then make a line segment bridging the break and combining all the line segments.

library(lwgeom)
library(tidyverse)
library(sf)

line1 <- lines %>% 
  st_sf() %>% 
  st_combine() %>% 
  st_sf() %>% 
  st_line_merge() %>% 
  st_cast("LINESTRING") 

### this is still some manual work which needs to be improved 
newline1 <- c(st_startpoint(line1[3,]), st_endpoint(line1[2,])) %>%
  st_combine() %>% 
  st_cast("MULTIPOINT") %>%
  st_union() %>%
  st_cast("LINESTRING") %>% 
  st_sf()

newline2 <- c(st_startpoint(line1[2,]), st_endpoint(line1[1,])) %>%
  st_combine() %>% 
  st_cast("MULTIPOINT") %>%
  st_union() %>%
  st_cast("LINESTRING") %>% 
  st_sf()

line1[nrow(line1)+1,]<-newline1
line1[nrow(line1)+1,]<-newline2
###

line1_uni <- line1 %>% 
  st_sf() %>% 
  st_combine() %>% 
  st_sf() %>% 
  st_line_merge() %>% 
  st_cast("LINESTRING") %>% 
  st_sf()

line1_uni

Simple feature collection with 1 feature and 0 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 2880 ymin: 340255.4 xmax: 2914.214 ymax: 340350.2
Projected CRS: MGI / Austria GK East
                        geometry
1 LINESTRING (2880 340255.4, ...

Connect Linestrings

2 Answers2