2

I want to select lines (bahnstrecken) within a polygon (bundesland). My goal is a selection shown in the following screenshot, which I did manually. Before I used two tools: I converted the polygons to lines and splitted the linies (bahnstrecken) with the new line-layer which was formerly a polygon. The result is a lines layer, which geoemtries end exactly at the polygon borders.

enter image description here

The problem is this: When I use different spatial query methods (intersects, contains, touches, fully contains) I never get the result I want. How can I reach the selection result shown in the first screenshot above? I think I should use the WITHIN method, but the results of the spatial selection tool appear quite strange to me:

INTERSECT:

enter image description here

TOUCHES:

enter image description here

OVERLAPS:

enter image description here

WITHIN:

enter image description here

1

I surmise from your comments and your layer names that you want to calculate the total length of railway lines in each State, so if you are open to a pyqgis solution I can suggest the following approach which does not require splitting or trying to create selections of your line features. This will simply calculate the parts of the line geometries which intersect each polygon region and sum their length.

You can run this this code directly with a polygon layer of States and a line layer of railway lines with out modifying or splitting either one.

A few notes:

  • It is important that both layers are in a projected CRS (which I see that both of yours already are).
  • I tested this with some data which looks similar to yours. A polygon layer with 17 regions and a line layer of roads with 3043 features. The code took around 25 seconds to complete so it's not instantaneous. I have used a spatial index to speed up the querying, but calculating the geometry intersection seems to be the computationally heavy part.

Paste the following code into an editor in the Python console and run. Make sure that the layer names match yours (I used layer names 'staaten' for states and 'bahnstrecken' for railways but if yours are different you will just need to change the first 2 lines of the script). Also please forgive my rusty German in the print statement at the end of the script.

line_layer = QgsProject().instance().mapLayersByName('bahnstrecken')[0]
polygon_layer = QgsProject().instance().mapLayersByName('staaten')[0]

line_feats = [f for f in line_layer.getFeatures()]

idx = QgsSpatialIndex(line_layer.getFeatures())

for p in polygon_layer.getFeatures():
    total_line_length = 0
    candidates = idx.intersects(p.geometry().boundingBox())
    candidate_feats = [f for f in line_feats if f.id() in candidates]
    for l in candidate_feats:
        intersect = l.geometry().intersection(p.geometry())
        total_line_length += intersect.length()
    length_km = round((total_line_length/1000), 2)
    print('Staat {}: enthält {}km Bahnstrecke'.format(p.id(), length_km))

My test data looks like this:

enter image description here

Example of output:

enter image description here

1
  • Thank you Ben! This works perfectly and the result makes sense. – Enzo Baldini Sep 22 '20 at 8:28
2

A feature has one geometry, for example a road feature has a line geometry. There are multi-geometries (e. g. a country that has some islands) but strictly speaking, they still count as one single geometry.

Spatial queries are queries, they ask if features have a special spatial relationship and depending on the true/false answer, add the features to the result set. For example for the predicate "within" it will return all features whose geometries lie completely inside the others.

You want not the whole lines but just specific parts of them. Since the query (or "selection") tools will look at the whole geometry of a feature and then return the whole feature, including its geometry, as it is, you are not able to cut away parts of the geometries using these tools.

Your goal is to actually change the geometries.

For that you can use the overlay tools. The operation you want here is called "clip". QGIS' documentation for the Vector overlay -> Clip tool says:

This algorithm clips a vector layer using the features of an additional polygon layer. Only the parts of the features in the Input layer that fall within the polygons of the Overlay layer will be added to the resulting layer.

8
  • I don't get you. Before I did the selection, I changed the geometry of the lines to prevent lines overlapping at the polygon borders. What I did was: 1: converting the polygons to lines 2: split the lines (bahnstrecken) with the region-lines. The result should be line geometries that don't overlap the polygon borders. You can see in the first screenshot that the lines reach only until the polygon border. So the WITHIN should give me all the lines within the polygon. But some are missing.... – Enzo Baldini Sep 18 '20 at 9:56
  • What I want to say is: The line features are already clipped as you suggest in your answer. – Enzo Baldini Sep 18 '20 at 10:47
  • Oh wow, I glanced over that bit of your explanation, sorry. My answer still will be the best way to do it. Spatial predicates are hard as you noticed. WITHIN will not include those that touch the border at any point (like where you cut them). INTERSECTS will include those that touch the border but go outside. Clip with the polygons! – bugmenot123 Sep 18 '20 at 10:56
  • 1
    I thought this way its easier to understand the problem. – Enzo Baldini Sep 18 '20 at 14:54
  • 1
    One approach would be to slightly buffer your border lines (like 1 meter or less) and then clip the network lines with that. Afterwards "within" should work as intended. I'll add that as possible answer later. – bugmenot123 Sep 19 '20 at 6:51

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.