Clip lines to grid using pyQgis

Question

I have created a regular grid of (square) polygons using Processing->Create grid. I would now like to iterate over my grid squares and clip an underlying road network within each grid square. I can do this from the GUI but was wondering how to do it using pyQgis (or console).

With l1 the layer of lines and l2 the grid My attempts have been along the lines of

iter=l2.getFeatures()
for f in iter:
    processing.runalg("qgis:clip", l1, f.attributes(), None)

But this fails with a wrong parameter... The third argument should be a polygon or vector. How do I convert?

Where am I going wrong?

Answer 1

With processing, 'qgis:clip' is not the adequate algorithm because it would produce one shapefile for each feature intersection. In this case, it's preferable to use 'saga:linepolygonintersection'. The complete code is:

import processing

l1 = processing.getObjectFromName('route')
l2 = processing.getObjectFromName('grid')

path = processing.runalg("saga:linepolygonintersection", l1, l2, 0, None)

clipped_line = QgsVectorLayer(path['INTERSECT'],
                              'clipped_line',
                              'ogr')

QgsMapLayerRegistry.instance().addMapLayer(clipped_line)

and after it was ran, at the Python Console of QGIS, with this shapefiles:

I got:

It can be observed that the algorithm probably has a bug because there are some not produced features.

To overcome this bad result it could be used this alternative code:

import processing

l1 = processing.getObjectFromName('route')
l2 = processing.getObjectFromName('grid')

l1_feats = [ feat for feat in l1.getFeatures() ]
l2_feats = [ feat for feat in l2.getFeatures() ]

geom_wkt = []

for feat_l1 in l1_feats:
    for feat_l2 in l2_feats:
        if feat_l1.geometry().intersects(feat_l2.geometry()):
            geom = feat_l1.geometry().intersection(feat_l2.geometry())
            geom_wkt.append(geom.exportToWkt())

epsg = l1.crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'clipped_line',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(geom_wkt)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(geom_wkt[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

This one produces the expected result:

Answer 2

You're almost there. You can make yours work by replacing the third parameter with a file containing a single polygon as described in the code below: If you find it's slow you can parallelise as described in Parallelising GIS operations in PyQGIS?

import processing
import os
input_file="/path/to/input_file.shp"
clip_polygons_file="/path/to/polygon_file.shp"
output_folder="/tmp/test/"
input_layer = QgsVectorLayer(input_file, "input file", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(input_layer)
l2  = QgsVectorLayer(clip_polygons_file, "clip_polys", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(l2)
l2_dp=input_layer.dataProvider()
EPSG_code=int(tile_layer_dp.crs().authid().split(":")[1])
tile_no=0
iter = l2.getFeatures()
for f in iter:
    print "Tile no: "+str(tile_no)
    tile_no+=1
    geom = f.geometry()
    clip_layer=QgsVectorLayer("Polygon?crs=epsg:"+str(EPSG_code)+\
    "&field=id:integer&index=yes","clip_polygon", "memory")
    clip_layer_dp = clip_layer.dataProvider()
    clip_layer.startEditing()
    clip_layer_feature = QgsFeature()
    clip_layer_feature.setGeometry(geom)
    (res, outFeats) = clip_layer_dp.addFeatures([clip_layer_feature])
    clip_layer.commitChanges()
    clip_file = os.path.join(output_folder,"tile_"+str(tile_no)+".shp")
    write_error = QgsVectorFileWriter.writeAsVectorFormat(clip_layer, \
    clip_file, "system", \
    QgsCoordinateReferenceSystem(EPSG_code), "ESRI Shapefile")
    QgsMapLayerRegistry.instance().addMapLayer(clip_layer)
    output_file = os.path.join(output_folder,str(tile_no)+".shp")
    processing.runalg("qgis:clip", input_file, clip_file, output_file)
    QgsMapLayerRegistry.instance().removeMapLayer(clip_layer.id())