2

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.

Improve this question
7
  • Related: gis.stackexchange.com/questions/168697/…
    – smiller
    Jun 6, 2019 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: docs.qgis.org/2.18/en/docs/user_manual/processing_algs/qgis/… docs.qgis.org/2.8/en/docs/user_manual/processing_algs/qgis/…
    – Ben W
    Jun 6, 2019 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 7, 2019 at 11:57
  • @smiler I check that before, but I could not adapt to my code.
    – tiago
    Jun 7, 2019 at 11:59
  • 2
    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, 2019 at 0:49

2 Answers 2

Reset to default
3

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.

example

Improve this answer
1
  • 1
    @M Bain, well done!
    – Taras
    Nov 7, 2019 at 6:07
1

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})
Improve this answer
1
  • Thank you @Benw, I'll try it.
    – tiago
    Jun 9, 2019 at 20:21

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.