Generating the skeleton of polygon using PyQGIS

Question

I develop a plugin for my project which includes some useful tools. I reviewed some answers in GIS SE about generating skeleton of polygon. I couldn't find any solution which used PyQGIS. HCMGIS plugin has a tool "Skeleton / Medail Axes" but it uses only processing tools.

I want to generate the skeleton of a polygon using pure PyQGIS without using Processing tools. After doing some digging I found voronoiDiagram method defined in QgsGeometry. I don't figure it out how to use it correctly.

I used voronoiDiagram in the following script, but it throws an exception: AttributeError: 'list' object has no attribute 'disjoint'.

What am I doing wrong? How can I improve/change this code to get skeleton of polygon?

My code is:

# get polygon layer
polygonLyr = QgsProject.instance().mapLayersByName("polygons")[0]

# create new layer
centerlineLyr =  QgsVectorLayer("Linestring?crs=EPSG:3857", "Centerlines", "memory")
QgsProject.instance().addMapLayer(centerlineLyr)

# layer has one polygon
polygonGeom = list(polygonLyr.getFeatures())[0].geometry()

# create voronoi lines
lineCollection = polygonGeom.voronoiDiagram(edgesOnly=True).asMultiPolyline()

#select lines and add them to layer
centerlineLyr.startEditing()
for line in lineCollection:
    if not line.disjoint(polygonGeom):
        f = QgsFeature()
        f.setGeometry(line)
        centerlineLyr.addFeature(f)
centerlineLyr.commitChanges()

Answer 1

I've added a script and some explanation about it.

1. You need polygon boundary for testing if voronoi line intersects boundary.
2. You need densified geometry. Because if polygon boundary has a long segment, skeleton will not be as desired (see the image). Densified geometry is just for generating voronoi lines. I used 10 meter. Depending on your data, you may need to change the value.

3. voronoiDiagram returns QgsGeometry. You have to convert it into geometry list for iterating over it. Since the result is a MultiLineString type, you need to use asGeometryCollection instead of asMultiPolyline. The reason of getting error is that asMultiPolyline returns list of list of QgsPointXY. But asGeometryCollection returns list of QgsGeometry (in this case list of LineGeometry)
4. If a line in geom_coll doesn't intersect polygon boundary and it is not outside polygon, that line is a part of skeleton.
5. After detecting lines of skeleton, convert separate lines into multiline. (You can add lines to layer as separate line features, if you want. I preferred multiline)

skeleton_layer =  QgsVectorLayer("Linestring?crs=EPSG:3857", "Skeleton", "memory")
polygon_layer = QgsProject.instance().mapLayersByName("polygons")[0]

for feature in polygon_layer.getFeatures():
    poly_geom = feature.geometry()    
    # (1)
    boundary = poly_geom.convertToType(1) # polygon boundary    
    # (2)
    dens_geom = poly_geom.densifyByDistance(10) # densified geometry     
    # (3) voronoi line collection 
    geom_coll = dens_geom.voronoiDiagram(edgesOnly=True).asGeometryCollection()    
    # (4) line collection for skeleton
    line_coll = [line.convertToType(1).asPolyline() for line in geom_coll \
            if not line.intersects(boundary) and not line.disjoint(poly_geom)]

    # (5) convert individual lines into multiline 
    m = QgsGeometry.fromMultiPolylineXY(line_coll)
    f = QgsFeature()
    f.setGeometry(m)
    skeleton_layer.dataProvider().addFeatures([f])

QgsProject.instance().addMapLayer(skeleton_layer)