Clipping point layer with polygon overlay using PyQGIS

Question

The following script works just fine and as intended:

import processing
import time

start_tot = time.time()

path = 'C:\\Users\\tlind\\Dropbox\\Documents\\'

ext = ".shp"

municipalities = ["bjuv", "bromölla", "burlöv", "båstad", "eslöv", "helsingborg", "hässleholm", "höganäs", "hörby", "hööör", "klippan", "kristianstad", "kävlinge", "landskrona", "lomma", "lund", "malmö", "osby", "perstorp", "simrishamn", "sjöbo", "skurup", "staffanstorp", "svalöv", "svedala", "tommelilla", "trelleborg", "vellinge", "ystad", "åstorp", "ängelholm", "örkelljunga", "östra göinge"] 

layer = "all_rs_400kv"

for value in municipalities:

    start = time.time()
    input_layer = path+layer+ext
    input_municipality = path+value+ext

    output_layer = path+layer+"_"+value+ext

    processing.run("native:clip", {'INPUT':input_layer, 'OVERLAY':input_municipality, 'DISSOLVE': False, 'OUTPUT':output_layer})

    iface.addVectorLayer(output_layer, '', 'ogr')

    end = time.time()
    print("Elapsed time for", value, ":", end-start, "seconds.")

end_tot = time.time()
print()
print("Total time elapsed:", (end_tot-start_tot)/60, "minutes.")

However, when adding a point layer as 'layer' it does not work. I'm getting the following error:

Traceback (most recent call last):
 File "C:\PROGRA~1\QGIS32~1.3\apps\Python39\lib\code.py", line 90, in runcode
   exec(code, self.locals)
  File "<input>", line 1, in (module)
  File "<string>", line 22, in (module)
  File "C:\PROGRA~1/QGIS32~1.3/apps/qgis/./python/plugins\processing\tools\general.py", line 108, in run
   return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:\PROGRA~1/QGIS32~1.3/apps/qgis/./python/plugins\processing\core\Processing.py", line 187, in runAlgorithm
   ret, results = execute(alg, parameters, context, feedback, catch_exceptions=False)
  File "C:\PROGRA~1/QGIS32~1.3/apps/qgis/./python/plugins\processing\gui\AlgorithmExecutor.py", line 72, in execute
   results, ok = alg.run(parameters, context, feedback, {}, False)
_core.QgsProcessingException: Could not write feature into OUTPUT

Any ideas?

Answer 1

Using a set of the following features

and the code below

import processing
import time

start_tot = time.time()

path = 'C:/Users/Taras/Downloads/data/'

ext = ".shp"

municipalities = ["municipality_1", "municipality_2"] 

layer = "points"

for value in municipalities:

    start = time.time()
    input_layer = path+layer+ext
    input_municipality = path+value+ext

    output_layer = path+layer+"_"+value+ext

    processing.run("native:clip", {'INPUT':input_layer, 'OVERLAY':input_municipality, 'DISSOLVE': False, 'OUTPUT':output_layer})

    iface.addVectorLayer(output_layer, '', 'ogr')

    end = time.time()
    print("Elapsed time for", value, ":", end-start, "seconds.")

end_tot = time.time()
print()
print("Total time elapsed:", (end_tot-start_tot)/60, "minutes.")

I could recreate the same error

Traceback (most recent call last):
  File "D:\Software\apps\Python39\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 22, in <module>
  File "D:\Software/apps/qgis/./python/plugins\processing\tools\general.py", line 108, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "D:\Software/apps/qgis/./python/plugins\processing\core\Processing.py", line 187, in runAlgorithm
    ret, results = execute(alg, parameters, context, feedback, catch_exceptions=False)
  File "D:\Software/apps/qgis/./python/plugins\processing\gui\AlgorithmExecutor.py", line 72, in execute
    results, ok = alg.run(parameters, context, feedback, {}, False)
_core.QgsProcessingException: Could not write feature into OUTPUT

The problem is about:

Feature creation error (OGR error: Attempt to write non-point (MULTIPOINT) geometry to point shapefile.)

The same error can be discovered when writing the temporal output to a shapefile:

alg_params = {
    'INPUT': input_layer,
    'OVERLAY': input_municipality,
    'OUTPUT': 'memory:'
    }
    
res = processing.run("native:clip", alg_params)['OUTPUT']

save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.driverName = "ESRI Shapefile"
save_options.fileEncoding = "UTF-8"
transform_context = QgsProject.instance().transformContext()
error = QgsVectorFileWriter.writeAsVectorFormatV3(res, output_layer, transform_context, save_options)

if error[0] == QgsVectorFileWriter.NoError:
    print("success!")
else:
    print(error)

that will lead to this error:

(7,
'Feature write errors:\nFeature creation error (OGR error: Attempt to write non-point (MULTIPOINT) geometry to point shapefile.)...
\nOnly 0 of 13 features written.', 'C:\\Users\\Taras\\Downloads\\data\\points_municipality_1.shp', 'points_municipality_1')

There are several ideas on how to overcome this error:

1. Using the GeoPackage as an output file format

   output_layer = join(path, layer + "_" + value + '.gpkg')
   

2. Applying additionally the "Multipart to singleparts" geoalgorithm

   alg_params = {
       'INPUT' : previous_output,
       'OUTPUT' : output_layer
       }
   processing.run("native:multiparttosingleparts", alg_params)  
   

There are also several relevant suggestions:

1. There is no need to have the 'DISSOLVE' parameter. Check the geoalgorithm with processing.algorithmHelp("native:clip"):

   Clip (native:clip)
   
   This algorithm clips a vector layer using the features of an additional polygon layer. Only the parts of the features in the Input layer that fall within the polygons of the Overlay layer will be added to the resulting layer.
   
   The attributes of the features are not modified, although properties such as the area or length of the features will be modified by the clipping operation. If such properties are stored as attributes, those attributes will have to be manually updated.
   
   ----------------
   Input parameters
   ----------------
   
   INPUT: Input layer
   
       Parameter type: QgsProcessingParameterFeatureSource
   
       Accepted data types:
           - str: layer ID
           - str: layer name
           - str: layer source
           - QgsProcessingFeatureSourceDefinition
           - QgsProperty
           - QgsVectorLayer
   
   OVERLAY: Overlay layer
   
       Parameter type: QgsProcessingParameterFeatureSource
   
       Accepted data types:
           - str: layer ID
           - str: layer name
           - str: layer source
           - QgsProcessingFeatureSourceDefinition
           - QgsProperty
           - QgsVectorLayer
   
   OUTPUT: Clipped
   
       Parameter type: QgsProcessingParameterFeatureSink
   
       Accepted data types:
           - str: destination vector file, e.g. 'd:/test.shp'
           - str: 'memory:' to store result in temporary memory layer
           - str: using vector provider ID prefix and destination URI, e.g. 'postgres:…' to store result in PostGIS table
           - QgsProcessingOutputLayerDefinition
           - QgsProperty
   
   ----------------
   Outputs
   ----------------
   
   OUTPUT:  <QgsProcessingOutputVectorLayer>
       Clipped
   

2. Applying the Python's native library os.path to handle file paths and names

3. Working with the try except block is also a common Python practice

   try:
       processing.run("native:clip", alg_params)
   except QgsProcessingException as err:
       print(err)
   

4. Last but not least make use of the spatial index for the input layer

So, after all, your code may look as follows:

import time
import processing
from os.path import join, normpath

start_tot = time.time()

path = normpath('C://Users//Taras//Downloads//data//')
layer = "points"
ext = ".shp"
input_layer = join(path, layer + ext)

municipalities = ["municipality_1", "municipality_2"] 

for value in municipalities:

    start = time.time()
    
    input_municipality = join(path, value + ext)

    alg_params1 = {
        'INPUT' : input_layer,
        'OVERLAY' : input_municipality,
        'OUTPUT' : 'memory:'
        }
    
    res1 = processing.run("native:clip", alg_params1)['OUTPUT']

    output_layer = join(path, layer + "_" + value + ext)
    
    alg_params2 = {
        'INPUT' : res1,
        'OUTPUT' : output_layer
        }
        
    res2 = processing.run("native:multiparttosingleparts", alg_params2)  
  
    iface.addVectorLayer(output_layer, '', 'ogr')
    
    end = time.time()
    print("Elapsed time for", value, ":", end-start, "seconds.")

end_tot = time.time()
print()
print("Total time elapsed:", (end_tot-start_tot)/60, "minutes.")

and get the output like this:

---

References:

• Point Sampling Tool Output Empty
• How to fix the error '_core.QgsProcessingException: There were errors executing the algorithm.'
• Quick workaround for PyQGIS error “file is not a directory” when saving files
• Export to shapefile from DXF throws an error: Attempt to write a non-point (MULTIPOINT) geometry to point shapefile