1

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?

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:

enter image description here

I got:

enter image description here

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:

enter image description here

1

You're almost there. You can make yours work by replacing the thrid parameter with a file containing a single polygon as described inthe 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())

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.