Python - gdal.Polygonize produce invalid geometries

Question

gdal.Polygonize is making invalid geometries and later in processing chain an error is reported stating that there are Ring Self-Intersections. Not one, quite a few so manual solution is not an option in my case. Is there any solution to this problem? I am using Python 2.7 and QGIS. Hot spot (ring self-intersection) is circled in red.

Question is the same like QGIS Raster to Polygons produce invalid geometries, but I am looking for some solution using Python and its libraries.

Answer 1

You may try using the v.clean GRASS tool from the Processing Toolbox.

Using the Python Console, the mandatory parameters for this tool are retrievable by directly calling the help:

import processing
processing.alghelp("grass7:v.clean")

which gives:

ALGORITHM: v.clean - Toolset for cleaning topology of vector map.
    input <ParameterVector>
    tool <ParameterSelection>
    threshold <ParameterNumber>
    GRASS_REGION_PARAMETER <ParameterExtent>
    GRASS_SNAP_TOLERANCE_PARAMETER <ParameterNumber>
    GRASS_MIN_AREA_PARAMETER <ParameterNumber>
    output <OutputVector>
    error <OutputVector>


tool(Cleaning tool)
    0 - break
    1 - snap
    2 - rmdangle
    3 - chdangle
    4 - rmbridge
    5 - chbridge
    6 - rmdupl
    7 - rmdac
    8 - bpol
    9 - prune
    10 - rmarea
    11 - rmline
    12 - rmsa

More in detail, I suggest using it with the "break" cleaning tool, which breaks lines at each intersection. Otherwise, other cleaning tools of interest could be (I also post them for the sake of clearness and mainly because I don't know your issue in detail):

• "snap": Snaps lines to vertex in threshold
• "rmdangle": Removes dangles
• "rmline": Removes all lines or boundaries of zero length, threshold is ignored
• "rmsa": Removes small angles between lines at nodes

Assuming that your polygonized layer is the active layer, you may directly run the algorithm from the Python Console:

import processing
layer = iface.activeLayer()

extent = layer.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

processing.runalg("grass7:v.clean", {"input": layer,
                    "tool": 0,
                    "threshold": 0.1,
                    "GRASS_REGION_PARAMETER":"%f,%f,%f,%f" %(xmin, xmax, ymin, ymax),
                    "GRASS_SNAP_TOLERANCE_PARAMETER": -1,
                    "GRASS_SNAP_TOLERANCE_PARAMETER": 0.0001,
                    "output": 'C:/Users/path_to_the_first_output/cleaned.shp',
                    "error": 'C:/Users/path_to_the_second_output/errors.shp'})

In this example, I set the parameter "tool" to 0 because I'm using the "break method": if you want to test the other cleaning tools, you just need to change the number with reference to the help of the algorithm I posted above.