1

I have a linestring layer that I want to clip with a polygon layer. So far so good. However, as these polygons are river catchments starting from a sampling site, they are sometimes overlapping (due to multiple sampling sites on a river). When this is the case, I want the clipped linestring to be duplicated according to the number of overlapping polygons (i.e. I want to have a linestring for each polygon within each polygon).

To illustrate you my example I've uploaded two sample shapes here. It contains 3 overlapping polygons and a linestring. My actual data consists of ~ 2500 polygons and a stream network for whole Germany.

enter image description here

One possibility that comes to my mind would be to loop through all the polygons' unique ID and clip them individually with the line network. However I ask myself, if there are already functions in QGIS, PostGIS, or R that can do this?

using the loop strategy in R I came with the code below. It's basically working, but I still have problems when a line has multiple segments within a polygon and it's extremely slow compared to the clip function in QGIS for example.

require(sf)
ezg_fin = st_read('polygons.shp')
gew_typ = st_read('line.shp')

todo = ezg_fin$site_id

line_list = list()
for (i in 1:length(todo)) {
  catchm = ezg_fin[ezg_fin$site_id == todo[i], ]
  line_sub = gew_typ[catchm, ] # 1st spatial subset the stream segment (fast!) - only possible 'cause stream network is segmented
  message(i, '. ', 'Clipping river at: ', todo[i])
  line_clip = st_intersection(line_sub, catchm) # 2nd clip only the subsetted segement (whole stream network would take way too long!)

  line_list[[i]] = line_clip
}

fin = do.call(rbind, line_list)
  • Did you already try anything? I would say postgis is very suitable for the job but it helps if you already have some code or workflow that we can comment on. – tilt Jun 1 '17 at 13:06
  • I'm currently trying it in R using a spatial subset [, as the actual line network is not unified (i.e. segments roughly the size of the bigger polygons). However this doesn't really clip the line at the border. I also tried sf::st_intersection() which is pretty slow 'cause my actual line network is really big. I elaborate a bit more and will provide some code then. My general strategy is looping through the polygon IDs. – andschar Jun 1 '17 at 13:12
  • You could use the sp classes, imported using rgdal::readOGR, and the sp::over function or intersect functions in rgeos. – Jeffrey Evans Jun 1 '17 at 16:47
4

There is a great wrapper function in the raster package, using rgeos::gIntersects, that will return exactly what you are after.

library(sp)
library(rgdal)
library(raster)

l <- rgdal::readOGR(getwd(), "line")
p <- rgdal::readOGR(getwd(), "polygons") 
sp::proj4string(p) <- sp::proj4string(l)

lp <- raster::intersect(l,p)
str(lp@data) 
  plot(p)
  plot(lp, col="red", add=TRUE)
2

QGIS example using Clip.

enter image description here

  1. Make sure your layers are on the same CRS. (I have reprojected polygons layer to EPSG:31467, same as line layer).
  2. Processing Toolbox | QGIS geoalgorithms | Vector overlay tools | Clip
  3. Select Input line and Clip polygons
  4. Click on Iterate over this layer button of the Clip layer
  5. Probably the output needs to be set to file ([Create temporary layer] may not work). Hit [Run]

One clipped line layer per polygon (total three) will be generated.

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.