TransWikia.com

How do I precisely select lines within a polygon in QGIS?

Geographic Information Systems Asked on December 16, 2020

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

2 Answers

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

Correct answer by Ben W on December 16, 2020

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.

Answered by bugmenot123 on December 16, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP