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.
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:
TOUCHES:
OVERLAPS:
WITHIN:
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:
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:
Example of output:
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
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP