Select by location using for-loop in Python Terminal

Question

I have two shapefiles ('drenagem' and 'rio_pri').

I need to loop through every feature from 'rio_pri' and save all those features from 'drenagem' that intersect 'rio_pri'.

So I wrote the following code to run in Python Terminal.

import processing
drenagem = QgsVectorLayer('/home/tiago/Documentos/Rios_mata_disol.shp','Rios_mata_disol','ogr')
rio_pri = QgsVectorLayer('/home/tiago/Documentos/rios_linha_costa.shp','rios_linha_costa','ogr')
# Loop
for r in rio_pri.getFeatures():
    processing.runalg("qgis:selectbylocation", drenagem, r,['intersects'],0,0)
    r_cod = r['COCURSODAG']
    path = r_cod+'.shp'
    processing.runalg('qgis:saveselectedfeatures', drenagem, path)

Then for every feature from 'rio_pri' I get this message at the console.

Unable to execute algorithm
Wrong parameter value: <qgis._core.QgsFeature object at 0x7ffb6620b770>

I guess I should not use the 'r' as an argument of qgis:selectbylocation, but I could not figure out what should I use instead. Any ideas?

I'm using QGIS 2.18.17 on Linux Mint 19.

Related: https://gis.stackexchange.com/questions/168697/performing-a-spatial-query-in-a-loop-in-pyqgis — SMiller, Jun 06 '19 at 20:02
Your main problem is that you are trying to run the processing algorithms for every feature in a layer! The algorithms themselves loop over the features. Just call select by location once, pass 'drenagem' as input and 'rio_pri' as intersect, then call save selected features on 'drenagem'. Check the docs: https://docs.qgis.org/2.18/en/docs/user_manual/processing_algs/qgis/vector_selection_tools.html#select-by-location https://docs.qgis.org/2.8/en/docs/user_manual/processing_algs/qgis/vector_general_tools/saveselectedfeatures.html — Ben W, Jun 06 '19 at 22:51
@BenW But I don't want to select features from 'drenagem' that intersect 'rio_pri' just ONCE. I need to do the select by location feature by feature from 'rio_pri'. — tiago, Jun 07 '19 at 11:57
@smiler I check that before, but I could not adapt to my code. — tiago, Jun 07 '19 at 11:59
OK, I misunderstood your goal. Do you mean that you want to save a new shapefile for every feature in rio_pri containing just the intersecting features from drenagem? What are the geometry types of the layers? — Ben W, Jun 07 '19 at 12:55
Hi @BenW, I want to build a shapefile containing all river that connect with a main river. To do that, starting from the main river selecting the first rivers that have intersections with the main one, save these selected. Then do another selection by location, starting from the last shapefile saved. Then repeat this proccess until all the river from the basin of the main river are together in one shapefile. Thanks for your help! — tiago, Jun 09 '19 at 20:21
Going by your last comment, I don't think my answer will help you much. An image or screenshot of what your existing data sets look like and what you want your output to look like would help. Making a sample of your data available (if possible) would probably be helpful as well. — Ben W, Jun 10 '19 at 00:49

score 3 · Answer 1 · edited Nov 07 '19 at 06:06

If your end goal is to create a separate shapefile for each river and all its tributaries. You could try this approach using a recursive query (@Taras check this out :) )

I am assuming the main rivers layer 'rios_linha_costa' has a column containing the name of the river called "river_name".

Create a Virtual Layer (through Layer > Add Layer > Add/Edit Virtual Layer...) with the query:

WITH RECURSIVE basin(river_name, stream_order, geometry) AS (
SELECT river_name, 1, geometry from rios_linha_costa
UNION 
SELECT river_name, stream_order + 1, Rios_mata_disol.geometry
FROM Rios_mata_disol, basin 
WHERE st_intersects(basin.geometry, Rios_mata_disol.geometry)
AND stream_order < 15
)

SELECT river_name, min(stream_order) AS stream_order, geometry
FROM basin
GROUP BY river_name, geometry

This will give you a new layer with two columns: "river_name" and "stream_order".

The line in the query AND stream_order < 15 is necessary to avoid the query infinitely looping by intersecting between adjacent and identical geometries - this could probably be made smarter. The performance was fine on my tiny test data. Feel free to adjust the value of 15 up or down to suit the depth of the stream hierarchy in your data.

You can then use the "Split vector layer" processing algorithm to save a separate file for each basin by specifying a river name as the unique identifier.

Ben W · Answer 2 · 2019-06-10T00:45:25.583

Edit: After reading your subsequent comments above about what you are trying to achieve, I now doubt this actually what you want...

I have written the code below which you should be able to adapt. I have tested this in QGIS 3.4 and have included some commented lines which should allow for adaptation to QGIS 2.18. No guarantees though as I have not used the QGIS 2 API for a while now. If it's possible for you, I would recommend upgrading to at least QGIS 3.4. The script below assumes the 2 layers are loaded in the current project. My output files are in geopackage format, named 'test_file_' appended with the id of the intersecting feature to which they belong- you can change as required. Also, double backslashes in the output file path string works on Windows- I'm not sure about Linux.

I hope this helps.

import processing

rio_pri = QgsProject().instance().mapLayersByName('rios_linha_costa')[0]
#QGIS 2# rio_pri = QgsMapLayerRegistry().instance().mapLayersByName('rios_linha_costa')[0]
drenagem = QgsProject().instance().mapLayersByName('Rios_mata_disol')[0]
#QGIS 2# drenagem = QgsMapLayerRegistry().instance().mapLayersByName('Rios_mata_disol')[0]
index = QgsSpatialIndex(drenagem.getFeatures())
for feat in rio_pri.getFeatures():
    select_ids = []
    geoms = feat.geometry()
    engine = QgsGeometry.createGeometryEngine(geoms.constGet())
    engine.prepareGeometry()
    intersection = index.intersects(geoms.boundingBox())
    candidates = drenagem.getFeatures(QgsFeatureRequest().setFilterFids(intersection))
    for f in candidates:
        if engine.intersects(f.geometry().constGet()):
            select_ids.append(f.id())
    drenagem.selectByIds(select_ids)
    out_layer = 'Path\\To\\Output\\Location\\test_file_{}.gpkg'.format(str(feat.id()))
    # change extension to .shp if preferred
    processing.run('native:saveselectedfeatures', {'INPUT': drenagem,
    'OUTPUT': out_layer})

Select by location using for-loop in Python Terminal

2 Answers2