12

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()
2
  • If you print(type(lineCollection)), it sends you back <class 'list'> – Vincent Bré Feb 21 '20 at 7:46
  • @VincentBré lineCollection is a list. I know that. The problem is items in the list. Each item of lineCollection is a list, too. I would expected that each item is a line. – menes Feb 21 '20 at 16:10
11
+50

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.

enter image description here

  1. 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)
  2. If a line in geom_coll doesn't intersect polygon boundary and it is not outside polygon, that line is a part of skeleton.
  3. 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)

enter image description here

4
  • Thank you. it works like a charm. My next step is to remove small dangles. Do you have any idea? – menes Feb 21 '20 at 16:18
  • You're welcome. I suggest you to search it on GIS SE, try something and ask another question if you stuck. – Kadir Şahbaz Feb 21 '20 at 16:31
  • i try to remove small dangles. still no complete success. – menes Dec 22 '20 at 1:10
  • I actually had tried when you asked. Unfortunately I couldn't find a proper solution to delete dangles at that time. I was stuck and I stopped trying :(. Removing dangles is a research subject all by itself. – Kadir Şahbaz Dec 22 '20 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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