3

I have a pyqgis script which customises the functionality of the grass algorithm r.lake.coords:

##Seed Points Lake=name
##DEM_input=raster
##Seed_vector_points=vector point
##Overspill_vector_boundary=vector polygon
##output=folder


# Import classes
from qgis.core import QgsRasterLayer
from qgis.utils import iface
from PyQt4.QtCore import QFileInfo
import processing

# Define layers
point_layer = processing.getObject(Seed_vector_points)
layer = processing.getObject(DEM_input)
boundary_layer = processing.getObject(Overspill_vector_boundary)
fileName = layer.source()
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(fileName, baseName)

# Get extent of the region
extent = iface.mapCanvas().extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

for feat in point_layer.getFeatures():
    geom = feat.geometry()
    attr = feat.attribute("Z_ADD_1")
    i = 0
    while i < 5: 
        output_1 = processing.runalg('grass:r.lake.coords', layer, "%f"%(attr), "%f,%f"% (geom.asPoint().x(), geom.asPoint().y()) , False, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), 0.00, None)
        output_2 = processing.runalg('saga:rastercalculator', output_1['lake'], '','a>0', False, 3, None)
        output_3 = processing.runalg('gdalogr:polygonize', output_2['RESULT'], "DN", None)
        output_4 = processing.runalg('qgis:dissolve',output_3['OUTPUT'], True, "DN", None)
        output_5 = processing.runalg("qgis:extractbylocation", output_4['OUTPUT'], boundary_layer, u'within', output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + "_%.2f"%(attr) + 'm.shp')
        i += 1
        attr += 0.5

My problem is that the final algorithm which the script calls upon produces additional blank shapefiles where there is no geometry, this is technically correct if the outputs from output_4 are not within the overspill boundary.

However, I would ideally like to have the script set up using some sort of python logic where it avoids producing final shapefiles with no geometry.

5

I prevent empty layers from loading with this code snippet:

layer = iface.activeLayer()
if not self.layer.isValid():
   iface.messageBar().pushMessage("Info:", "Layer not valide", level=QgsMessageBar.WARNING)
   pass
elif layer.featureCount() == 0:
   pass 
elif layer.featureCount() != 0:
   print "do something usefull here"
   QgsMapLayerRegistry.instance().addMapLayer(layer)

In your use case, you would iterate over the 5 outputs and and check if the layer is empty and then load the layers with QgsMapLayerRegistry.instance().addMapLayer(layer) into QGIS

1
  • Great code for preventing the loading of null geometries into QGIS @LaughU, I'll certainly be using this in the future :)
    – Ian
    Aug 31 '17 at 13:59
3

Great answer by @LaughU! If you want to prevent processed layers with no geometry to have no output, you could use the None parameter to save the output as a temporary file initially, check if this has any features and if so, write to a shapefile otherwise skip it:

...
output_5 = processing.runalg("qgis:extractbylocation", output_4['OUTPUT'], boundary_layer, u'within', 0.00, None)
layer = processing.getObject(output_5['OUTPUT'])
if not layer.featureCount() == 0:
    QgsVectorFileWriter.writeAsVectorFormat(layer, output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + "_%.2f"%(attr) + 'm.shp', "utf-8", None, "ESRI Shapefile")
    i += 1
    attr += 0.5
3
  • Hi Joseph, it seems to be working fine for a single iteration (I assume) as I get one output and then I get the error message "'NoneType' object has no attribute 'getitem' See log for more details". Having tested this out manually, I should be getting two valid results.
    – Ian
    Aug 31 '17 at 14:01
  • In the console the same error is "Error: Wrong parameter value: <qgis._core.QgsVectorLayer object at 0x155228A0>"
    – Ian
    Aug 31 '17 at 14:06
  • @Ian - Apologies, I missed the precision parameter in the algorithm. This should avoid the wrong parameter error atleast. Also, are you sure i = 0 should not be outside the loop? Because it is resetting to 0 for every feature.
    – Joseph
    Aug 31 '17 at 14:14
1

Having used both of your excellent answers, I have come up with the following solution which works perfectly:

##Seed Points Lake=name
##DEM_input=raster
##Seed_vector_points=vector point
##Overspill_vector_boundary=vector polygon
##output=folder


# Import classes
from qgis.core import (QgsRasterLayer, QgsVectorFileWriter)
from qgis.utils import iface
from qgis.gui import QgsMessageBar
from PyQt4.QtCore import QFileInfo
import processing

# Define layers
point_layer = processing.getObject(Seed_vector_points)
layer = processing.getObject(DEM_input)
boundary_layer = processing.getObject(Overspill_vector_boundary)
fileName = layer.source()
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(fileName, baseName)


# Get extent of the region
extent = iface.mapCanvas().extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

for feat in point_layer.getFeatures():
    geom = feat.geometry()
    attr = feat.attribute("Z_ADD_1")
    i = 0
    while i < 5: 
        output_1 = processing.runalg('grass:r.lake.coords', layer, "%f"%(attr), "%f,%f"% (geom.asPoint().x(), geom.asPoint().y()) , False, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), 0.00, None)
        output_2 = processing.runalg('saga:rastercalculator', output_1['lake'], '','a>0', False, 3, None)
        output_3 = processing.runalg('gdalogr:polygonize', output_2['RESULT'], "DN", None)
        output_4 = processing.runalg('qgis:dissolve',output_3['OUTPUT'], True, "DN", None)
        output_5 = processing.runalg("qgis:extractbylocation", output_4['OUTPUT'], boundary_layer, u'within', None)
        layer_1 = processing.getObject(output_5['OUTPUT'])
        if not layer_1.featureCount() == 0:
            QgsVectorFileWriter.writeAsVectorFormat(layer_1, output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + "_%.2f"%(attr) + 'm.shp', "utf-8", None, "ESRI Shapefile")
        else:
            iface.messageBar().pushMessage("Info:", "Process Complete :)", level=QgsMessageBar.SUCCESS)
        i += 1
        attr += 0.5

Thanks again for the solutions guys, much appreciated.

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.