"Lines to Polygon: output missing features when saving to SHP and input is multipart

[Muesgen]

What is the bug or the crash?

I have found a bug with the "lines to polygon" tool of QGIS and would like to know if the bug is known and if there is a workaround?

I want to convert my MultiLineString Z to MultiPolygon Z with Qgis inside a plugin. The results from the processor tool are different from the Python console than from the tool interface.

Polygons are missing in the layer created via the Python console, while all polygons are in the output via the interface. However, if I do not select a temporary output but a layer file to be written out, the polygon geometry is also missing via the interface .

Steps to reproduce the issue

See

Versions

QGIS 3.18 and Python 3.7

Supported QGIS version

• I'm running a supported QGIS version according to the roadmap.

New profile

• I tried with a new QGIS profile

Additional context

No response

[gioman]

QGIS 3.18 and Python 3.7

@Muesgen please test on 3.22.3, 3.18 is unsupported. Also add here all the details you have added on SE (where someone has left an answer with a clue of that could be the problem).

[Muesgen]

I tested on 3.22.3 Version and the problem is still there.

You can access the example data from my Google Drive in shapefile format: https://drive.google.com/drive/folders/1I_uhDnlvqjhrHZfUWmIpZoi-C5nDTG47?usp=sharing

The lines are extracted from citygml file from this source: https://www.opengeodata.nrw.de/produkte/geobasis/3dg/lod2_gml/lod2_gml/LoD2_32_352_5662_1_NW.gml

Python Console executed inside QGIS. 3.22.3

test = "c:/Users/abc/Desktop/lod2_dev_test/LoD2_32_352_5662_1_NW.shp"
out = "c:/Users/abc/Desktop/lod2_dev_test/Polygone.shp"
processing.run("qgis:linestopolygons", {'INPUT':test,'OUTPUT':out}) 

The output looks like this:

if I run the same algorithm with the interface, the result is correct:

Here the recorded command from the protocol:
processing.run("qgis:linestopolygons", {'INPUT':'c:/Users/abc/Desktop/lod2_dev_test/LoD2_32_352_5662_1_NW.shp','OUTPUT':'TEMPORARY_OUTPUT'})

The Polygons look like this:

The lines definitely meet, that I checked in the WKT format of the shape.

On StackExchange there was an answer that the error happens when writing out the shapefile, because with GeoPackage there was no problem with writing out the correct geometry.

Quote frome SE:

It seems to be an issue with using shapefile as output. When using GeoPackage, it works seamlessly.

test = "c:/Users/abc/Desktop/LoD2_32_352_5662_1_NW.shp"
out = "ogr:dbname='C:/Users/abc/Desktop/out.gpkg' table='Polygone' (geom)"

# if you use the next line, you loose a geometry
# out = 'C:/Users/abc/Desktop/Polygone.shp'

result = processing.run("qgis:linestopolygons",
                        {'INPUT':test,'OUTPUT':out})["OUTPUT"]

[gioman]

@Muesgen

1. there is no difference on how the tool runs from the pyqgis console or the GUI. The difference you observe is because with the GUI the default output format is GPKG. If you force the output to SHP ypu get the same result

2. the problem seems to have origin from the fact that your lines are 1 unique multi-geometry. if you explode it to singleparts then you get the right result even when exporting to SHP

[agiudiceandrea]

@gioman it

2. the problem seems to have origin from the fact that your lines are 1 unique multi-geometry. if you explode it to singleparts then you get the right result even when exporting to SHP

@gioman, while I agree it is unusual that the layer contains only 1 multipart feature geometry (MultiLineStringZ) with 9871 parts (most of them are also duplicates: only 4811 of them are unique geometries), it seems to me this is not the real cause of the issue.

First of all, the "Lines to polygons" (qgis:linestopolygons) algorithm has code in place to supposedly correctly handle multipart (QgsGeometryCollection) geometries:

QGIS/python/plugins/processing/algs/qgis/LinesToPolygons.py

Lines 123 to 136 in bbfa18f

	def getSurfaces(self, geometry):
	surfaces = []
	if isinstance(geometry, QgsGeometryCollection):
	# collection
	for i in range(geometry.numGeometries()):
	surfaces.extend(self.getSurfaces(geometry.geometryN(i)))
	else:
	# not collection
	if geometry.vertexCount() > 2:
	surface = QgsPolygon()
	surface.setExteriorRing(geometry.clone())
	surfaces.append(surface)
	return surfaces

Second, there is an inconsistency: exporting to a Shapefile layer produces different results than exporting to a GeoPackage or a Memory layer.

Moreover, as you have already find, if you convert ("Multipart to singleparts") the original layer to a layer with 9871 singleparts features, then you get the same results applying the "Lines to polygons" alg to this layer and exporting to both GeoPackage and Shapefile.
However, if you apply the "Collect geometries" alg to the resulting layer (to have a multipart layer) and export it to a Shapefile layer, then the exported Shapefile layer will have again "missing" polygons, while there will be no missing polygons when exporting to a GeoPackage layer.

Is it possible to replicate the issue using a subsets of the original data, which contain 1 multipart feature (MultiLineStringZ) with only 7 parts: Issue47288_Lines.zip

Applying the "Lines to polygons" alg and exporting to a GeoPackage or a Memory layer (Issue47288_Polygons.gpkg in Issue47288_Polygons_GPKG.zip):

Applying the "Lines to polygons" alg and exporting to a Shapefile layer (Issue47288_Polygons.shp in Issue47288_Polygons_SHP.zip):

ogrinfo Issue47288_Lines.shp Issue47288_Lines -geom=SUMMARY reports:

Geometry: 3D Line String
Feature Count: 1
...
OGRFeature(Issue47288_Lines):0
  FID (Integer64) = 0
  MULTILINESTRING : 7 geometries:
LINESTRING : 16 points
LINESTRING : 5 points
LINESTRING : 5 points
LINESTRING : 5 points
LINESTRING : 10 points
LINESTRING : 5 points
LINESTRING : 5 points

ogrinfo Issue47288_Polygons.gpkg Issue47288_Polygons -geom=SUMMARY reports:

Geometry: 3D Multi Polygon
Feature Count: 1
OGRFeature(Issue47288_Polygons):0
  MULTIPOLYGON : 7 geometries:
POLYGON : 16 points
POLYGON : 5 points
POLYGON : 5 points
POLYGON : 5 points
POLYGON : 10 points
POLYGON : 5 points
POLYGON : 5 points

ogrinfo Issue47288_Polygons.shp Issue47288_Polygons -geom=SUMMARY reports:

Geometry: 3D Polygon
Feature Count: 1
...
OGRFeature(Issue47288_Polygons):0
  FID (Integer64) = 0
  MULTIPOLYGON : 2 geometries:
POLYGON : 16 points, 2 inner rings (5 points, 5 points)
POLYGON : 10 points, 3 inner rings (5 points, 5 points, 5 points)

I think the issue is related to the GDAL/OGR library which QGIS relies on to write and read Shapefile format layers or some limitation of the Shapefile format.

Tested with QGIS 3.16.16, 3.22.3, 3.23.0 from OSGeo4W v2 (GDAL 3.4.1).

[gioman]

I think the issue is related to the GDAL/OGR library which QGIS relies on to write and read Shapefile format layers or some limitation of the Shapefile format.

@agiudiceandrea thanks for the analysis! I'm fine with the conclusions :)

[agiudiceandrea]

This may be eventually fixed by OSGeo/gdal#7157.

[agiudiceandrea]

The issue has been fixed in the GDAL/OGR library with PR OSGeo/gdal#7157: it will not occur using a GDAL/OGR release version >= 3.7.0. On Windows, official OSGeo4W stand-alone and network installers of QGIS 3.28.7 and QGIS 3.30.3 and later versions will use GDAL/OGR release version >= 3.7.0.

[alexbruy]

Seems can be closed now as upstream issue was fixed.

[agiudiceandrea]

While the issue no longer occurs using the latest QGIS versions for Windows distributed by OSGeo4W (which currently ships GDAL 3.7.2), the issue still occurs e.g. using the latest QGIS versions for Ubuntu Lunar and Jammy (GDAL 3.4.1 and 3.6.2) and for Debian bookworm (GDAL 3.6.2).

[alexbruy]

I'm aware of that. But upstream issue is fixed and we don't need to do anything on the QGIS side. How keeping this open will help users with old GDAL?