1

I'm developing a standalone script in PyQGIS and have downloaded data from OpenStreetMaps via Overpass Turbo. The data is given as GeoJSON and has been converted to gpkg. It contains multiple geometries (LineString, Polygon, Points). As I need to access each geometry individually, I want each geometry type saved as a single layer. The code is running without any issues, however, no output is produced. I am new to Python.

inputDir = "C:\\Users"
outDir = "C:\\Users\\temp"

def separate_by_geometry_type(inputDir, outDir):
    input_layer = QgsVectorLayer(inputDir, 'input_layer', 'ogr')

    # Group the features by geometry type
    groups = {}
    for feature in input_layer.getFeatures():
        geom_type = feature.geometry().type()
        if geom_type not in groups:
            groups[geom_type] = []
        groups[geom_type].append(feature)

    # Save each group to a separate layer
    for geom_type, features in groups.items():
        # Create a new memory layer
        memory_layer = QgsVectorLayer(f'{geom_type}?crs={input_layer.crs().authid()}', f'{str(geom_type).lower()}_layer', 'memory')

        # Add features to the memory layer
        memory_layer_data = memory_layer.dataProvider()
        memory_layer_data.addFeatures(features)

        # Save the memory layer to a GeoPackage
        output_layer_path = f'{outDir}/{str(geom_type).lower()}_layer.gpkg'
        coordinateTransformContext=QgsProject.instance().transformContext()
        options = QgsVectorFileWriter.SaveVectorOptions()
        options.driverName = 'GPKG'
        QgsVectorFileWriter.writeAsVectorFormatV3(memory_layer, output_layer_path, coordinateTransformContext, options)

separate_by_geometry_type(inputDir, outDir)
2
  • Either way is fine for me. My goal is that python automatically detects whether a layer has multiple geometries and saves the polygons, lines and points in an individual layer. If I understand it correctly, the tool splits multipart into singlepart features, I, however, want each geometry type layer to remain its multipart feature property.
    – user210818
    Commented Dec 12, 2023 at 8:47
  • Sure. Here is a link to a layer I used. Thank you!
    – user210818
    Commented Dec 12, 2023 at 10:45

1 Answer 1

2

There are probably plenty of approaches available to tackle your issue, especially based on pure Python and preparation of the OpenStreetMap data. However, let's stick with PyQGIS and its troubleshooting.

The most relevant tool in scope of PyQGIS is the "Filter by geometry type", that filters features by their geometry type and provides outputs based on whether they have a point, line or polygon geometry.

Even, if you adjust your code correspondingly:

from os.path import abspath, join
from qgis import processing

inputFile = abspath("D:/qgis_test/OSM_archaeological_sites.geojson")
outputDir = abspath("D:/qgis_test/temp")

params = {
        'INPUT' : inputFile,
        'POINTS' : join(outputDir, "point_layer.gpkg"),
        'LINES' : join(outputDir, "polyline_layer.gpkg"),
        'POLYGONS' : join(outputDir, "polygon_layer.gpkg"),
        'NO_GEOMETRY' : 'TEMPORARY_OUTPUT'
        }

processing.run("native:filterbygeometry", params)

you may encounter the following exception:

Traceback (most recent call last):
  File "C:\OSGeo4W\apps\Python39\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 19, in <module>
  File "<string>", line 14, in separate_by_geometry_type
  File "C:\OSGeo4W/apps/qgis-ltr/./python/plugins\processing\tools\general.py", line 108, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:\OSGeo4W/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 187, in runAlgorithm
    ret, results = execute(alg, parameters, context, feedback, catch_exceptions=False)
  File "C:\OSGeo4W/apps/qgis-ltr/./python/plugins\processing\gui\AlgorithmExecutor.py", line 72, in execute
    results, ok = alg.run(parameters, context, feedback, {}, False)
_core.QgsProcessingException: Could not create layer D:\qgis_test\temp\point_layer.gpkg: Unsupported type for field relations

As you can see it is happening because the input GeoPackage file contains a certain field "@relations" (something vital for OpenStreetMap data), which causes the breakdown. For example:

"@relations": [
          {
            "role": "building",
            "rel": 122481,
            "reltags": {
              "historic": "archaeological_site",
              "name": "Eppinger Linie",
              "ref": "Eppinger Linien",
              "type": "collection"
            }
          }
        ]

It is a list with a dict inside. One can also notice it after uploading the original file into QGIS, which, by default, interprets it as a Key/Value-field. So, the issue is not about export, it is about how QGIS treats certain fields while importing/reading them.

Please, keep in mind, that this code may affect your original file.

# imports
from os.path import abspath, join
from qgis.core import edit, QgsVectorLayer, QgsField
from PyQt5.QtCore import QVariant
from qgis import processing

inputFile = abspath("D:/qgis_test/OSM_archaeological_sites.geojson")
outputDir = abspath("D:/qgis_test/temp")

layer = QgsVectorLayer(inputFile, "mygeojson", "ogr")
   
with edit(layer):
    layer.addAttribute(QgsField("relations", QVariant.String))
    layer.updateFields()

    issue_field = "@relations"

    for feat in layer.getFeatures():
        feat["relations"] = str(feat[issue_field])
        layer.updateFeature(feat)

    indx = layer.fields().indexFromName(issue_field)
    layer.deleteAttribute(indx)

params = {
        'INPUT' : layer,
        'POINTS' : join(outputDir, "point_layer.parquet"),
        'LINES' : join(outputDir, "polyline_layer.parquet"),
        'POLYGONS' : join(outputDir, "polygon_layer.parquet"),
        'NO_GEOMETRY' : 'TEMPORARY_OUTPUT'
        }

processing.run("native:filterbygeometry", params)

As you can see, the above code will produce a feasible output in GeoParquet file format. Also works for GeoJSON.

result1

However, this code will not work for GeoPackage, it yields the following exception:

Traceback (most recent call last):
  File "C:\OSGeo4W\apps\Python39\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 31, in <module>
  File "C:\OSGeo4W/apps/qgis-ltr/./python/plugins\processing\tools\general.py", line 108, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:\OSGeo4W/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 187, in runAlgorithm
    ret, results = execute(alg, parameters, context, feedback, catch_exceptions=False)
  File "C:\OSGeo4W/apps/qgis-ltr/./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 POLYGONS

The same for POINTS, and LINES. Perhaps someone can explain better why it happens. I suspect that GeoPackage does not accept a MultiPolygon to be written as a Polygon.

With a hint from @Pieter, a promising attempt for GeoPackage can be made like this:

# imports
from os.path import abspath, join
from qgis.core import edit, QgsVectorLayer, QgsField, QgsProject, QgsVectorFileWriter
from PyQt5.QtCore import QVariant
from qgis import processing
from qgis.utils import iface

inputFile = abspath("D:/qgis_test/OSM_archaeological_sites.geojson")
outputDir = abspath("D:/qgis_test/temp")

layer = QgsVectorLayer(inputFile, "originalgeojson", "ogr")

with edit(layer):
    layer.addAttribute(QgsField("relations", QVariant.String))
    layer.updateFields()

    issue_field = "@relations"

    for feat in layer.getFeatures():
        feat["relations"] = str(feat[issue_field])
        layer.updateFeature(feat)

    indx = layer.fields().indexFromName(issue_field)
    layer.deleteAttribute(indx)

params = {
    'INPUT': layer,
    'POINTS': 'TEMPORARY_OUTPUT',
    'LINES': 'TEMPORARY_OUTPUT',
    'POLYGONS': 'TEMPORARY_OUTPUT',
    'NO_GEOMETRY': 'TEMPORARY_OUTPUT'
}

result = processing.run("native:filterbygeometry", params)

for key, value in result.items():
    if isinstance(value, QgsVectorLayer) and key != 'NO_GEOMETRY':
        layer1 = result.get(key)

        params1 = {
            'INPUT': layer1,
            'OUTPUT': 'TEMPORARY_OUTPUT'
        }

        layer2 = processing.run("native:promotetomulti", params1)['OUTPUT']

        params2 = {
            'INPUT': layer1,
            'FIELD_NAME': "fid",
            'START': 1,
            'GROUP_FIELDS': [],
            'MODULUS': 0,
            'SORT_ASCENDING': True,
            'SORT_EXPRESSION': '',
            'SORT_NULLS_FIRST': False,
            'OUTPUT': 'memory:'
        }

        layer3 = processing.run("native:addautoincrementalfield", params2)['OUTPUT']

        if layer3.isValid():
            output_layer_path = join(outputDir, f"layer_{layer.wkbType()}.gpkg")

            coordinateTransformContext = QgsProject.instance().transformContext()
            save_options = QgsVectorFileWriter.SaveVectorOptions()
            save_options.driverName = "GPKG"
            save_options.fileEncoding = "UTF-8"
            error = QgsVectorFileWriter.writeAsVectorFormatV3(
                layer=layer3,
                fileName=output_layer_path,
                transformContext=coordinateTransformContext,
                options=save_options)

            if error[0] == QgsVectorFileWriter.NoError:
                print("Layer successfully created and saved!")
                iface.addVectorLayer(output_layer_path, 'try', 'ogr')
            else:
                print(error[1])

However, it leads to this error:

Feature write errors:
Feature creation error (OGR error: sqlite3_exec(CREATE TABLE "layer_6" (...) failed: duplicate column name: fixme)

because there are two fields with the same name: "FIXME" and "fixme". It should be either merged, removed or renamed. They contain the same information.

For the moment I have deleted "FIXME"-field:

indx2 = layer.fields().indexFromName("FIXME")
layer.deleteAttribute(indx2)

After the above code should work properly and produce such output:

result2

FIN


References:

2
  • Geopackage indeed doesn't support mixed "single" and "multi" features in one layer. The typical solution to this is to convert all features to "multi". So I would try to use the "Promote to multipart" tool before doing "filterbygeometry".
    – Pieter
    Commented Dec 15, 2023 at 14:57
  • Also not super elegant, but an alternative that is easily scriptable could be to first do the "filterbygeometry" to geojson files, and then convert those using "native:promotetomulti" to gpkg?
    – Pieter
    Commented Dec 15, 2023 at 15:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.