I have a shapefile of river segments that I've imported into R. The segments range in length from 5-115km. I want to split the segments that are over 10km into 10km chunks. So I use st_segmentize. But this doesn't actually split a long line into multiple new shorter lines, which is what I want... Is there a way to do this? I think st_segmentize just adds vertices every 10km. So I just need to split the lines at the vertices.

Code so far:


rivers = st_read("./YukonMain4_copy.shp", layer = 'YukonMain4_copy')

seg = st_segmentize(rivers, units::set_units(10,km))

Re. my above comment

The issue is, my resulting lines have lengths mostly from 200-500m, and then there are random onesthat are 700000m long. Plus, this data set is of multiple rivers, so now I have lines connecting places that should not be connected. Any ideas?

Here is my code:

  1. Import file, segmentize the data, and cast it into the point cloud
  2. Create an index to combine segments by (hopefully so that they are approx 10km in length
  3. Combine things back into lines.

    rivers=st_read("./YukonMain4_copy.shp", layer='YukonMain4_copy')
    seg = st_segmentize(rivers, units::set_units(10000,m))
    seg2= st_cast(seg, 'POINT')
    length = length(seg2$OBJECTID)/2
    index = rep(1:length, each=2)
    rivers.index=cbind(seg2, index)
    to_line <- function(points) st_cast(st_combine(points), "LINESTRING") %>% .[[1]] 
    rivers.nest = rivers.index %>% group_by(index) %>% nest
    segments = rivers.nest %>% pull(data) %>% map(to_line)%>%st_sfc(crs = 4326)
    lines = rivers.nest %>% select(data)%>% st_sf(geometry=segments)
  • Have you tried simply defining the dfMaxLength argument in the units of the projection (eg., dfMaxLength=10000)? Dec 8 '17 at 17:47
  • I have, but it doesn't seem to do anything. My rivers variable has 87 linesegments, and the seg variable still has 87 segments, and if I use st_length to check the lengths of the all the segments, they are identical.
    – Ana
    Dec 8 '17 at 17:53
  • The epsg (coordinate system) is 4326, which is the WGS84. So I think my data is in geographic coordinates. Oddly enough, when I use `st_length' I get units of meters... Do I need to add a projected coordinate system?
    – Ana
    Dec 8 '17 at 17:56
  • 1
    st_segmentize turns LINESTRINGs into LINESTRINGs, but they have more vertexes. If you plot(st_cast(segmentedthing, "POINT")) you should see those extra points. Next step is, I guess, to rebuild small segments from pairs of those points with the right attributes, but at this point (pun not intended) I wonder if there's another way...
    – Spacedman
    Dec 8 '17 at 17:56
  • 1
    See "Variations along a storm track" in this blog post. This may help to get the individual line segments extracted.
    – TimSalabim
    Dec 8 '17 at 18:05

I was looking for exactly this function to use with NHDPlus data. I spent some time with JMT2080AD's functions but needed a bit more scalable solution.

The function in this gist does the trick pretty nicely. There is a reproducible example at the bottom of the gist too.

--- Edit:

There is also now a function to do this in the lwgeom package. https://github.com/r-spatial/lwgeom/issues/16


I have a function I made for this. It uses sp, but it might be what you need. Use projected coordinates.



Well, I found some bugs in my function, so that's a bummer. It'll take me some time to sort that out. When I applied it to a large data set of numerous rivers, I found that r's seq function has some floating point issues and isn't working as I expected. That said, this works on roughly 95% of the lines I pass to it. If you have a study area in mind this might be ok, as it as been for me in the past.

This example assumes that you can work with columns in data.table that have spatial objects. There is probably a better solution out there than this. I haven't had to apply this to such a big file before.

## shapefile -> https://pubs.usgs.gov/dds/dds-81/TopographicData/River/
## river.shp, not UTM. Using geographic coordinates to match OP's example.


## reading in river.shp using sf and converting to UTM
riv <- st_read("./rivers.shp")
riv <- st_transform(riv, 3157) 

## converting to sp and making a list of river segements by record in spatial df
rivSegs <- as(riv, "Spatial")
rivSegs <- lapply(1:length(rivSegs), function(x) rivSegs[x,])

## setting up parallel cluster for iterating over spatial object
## this should work in windows or linux (mac also, probably..)
cl <- makeCluster(detectCores() - 1, outfile = "")
clusterExport(cl, ls())
clusterEvalQ(cl, {source("./polylineSplitter.r")})

## running splitLines against river segment list 
rivSplit <- parLapply(cl, rivSegs[1:40], splitLines, dist = 100)

## test plot
plot(rivSegs[[1]], add = T, col = 'red')

## building data.table of original attributes and adding sp list items
## converting back to sf for consistancy
out <- data.table(riv[1:40,])
out[,segments:=lapply(rivSplit, as, "sf")]
  • Apologies for the delayed reply. That looks helpful, thank you! I'm working with a global data set so projected coordinates is challenging but workable! It seems like the function works well on one line at a time. So I just have to loop through each line in the SpatialLinesDataFrame to split one line at a time. However, I would like to combine all the split lines in to a new SpatialLinesdataframe in the loop if that makes sense. But I can't seem to use like rbind or cbind to do this. Any tips?
    – Ana
    Dec 14 '17 at 16:00
  • Yes, that can be challenging but do (spdf behavior is weird using rbind). I will edit my post later today with an example.
    – JMT2080AD
    Dec 14 '17 at 17:15
  • Apologies for the delayed reply, but I just wanted to say thank you for the detailed example. That works well.
    – Ana
    Dec 18 '17 at 19:30
  • Oh great, I'm glad it works! Like I said, I've never tried it against so many records at once, just one segment at a time. This was a good test of that function. I see that I still have some work to do on it.
    – JMT2080AD
    Dec 18 '17 at 23:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.