Index de l'article

Centerlines (Skeleton)

This involves extracting the central lines of polygons (the skeleton). Here we use the plugins Centerlines from Geometric Attributes.

But it is not easy to get something satisfactory, even less when the polygons are of complex geometries and/or large sizes.

Here is a proposal where we create before a Clip of the polygons on a Grid before extracting the skeleton.

But you will have to play on the Grid and Centerlines parameters according to your context.

import shutil
import time
 
myPath = r'C:\\Users\\Georges\\Downloads\\_Test\\'
 
project = QgsProject.instance()
project.removeAllMapLayers()
project.clear()
iface.mapCanvas().refresh()
 
project = QgsProject.instance()
project.removeAllMapLayers()
project.clear()
iface.mapCanvas().refresh()
 
test = QgsVectorLayer(myPath + 'test.shp', 'test', 'ogr')
project.addMapLayer(test)
myTest = project.mapLayersByName('test')[0]
 
iface.mapCanvas().refresh()
 
def myFunc():
 
    myTest.dataProvider().createSpatialIndex()
 
    # PROCESS ÉTENDUE
    _test_etendue = myPath + 'test_etendue'
 
    if os.path.isdir(_test_etendue) == True:
        shutil.rmtree(_test_etendue)
    if os.path.isdir(_test_etendue) == False:
        os.mkdir(_test_etendue)
 
    test_etendue_path = _test_etendue + r'\\test_etendue.shp'
 
    processing.run("native:polygonfromlayerextent", {
        'INPUT':test,
        'ROUND_TO':0,
        'OUTPUT':test_etendue_path
    })
 
    test_etendue = QgsVectorLayer(test_etendue_path, 'test_etendue', "ogr")
    project.addMapLayer(test_etendue)
 
    myTestEtendue = project.mapLayersByName('test_etendue')[0]
    myTestEtendue.dataProvider().createSpatialIndex()
 
    # RÉCUPÉRER ÉTENDUE
    extent_test_etendue = test_etendue.extent()
    proj_test_etendue = test_etendue.crs()
    print(extent_test_etendue)
    print(proj_test_etendue)
 
    # PROCESS GRILLE
    _test_grille = myPath + 'grille'
 
    if os.path.isdir(_test_grille) == True:
        shutil.rmtree(_test_grille)
    if os.path.isdir(_test_grille) == False:
        os.mkdir(_test_grille)
 
    test_grille_path = _test_grille + r'\\test_grille.shp'
 
    processing.run("native:creategrid", {
        'TYPE':2,
        'EXTENT': extent_test_etendue,
        'HSPACING':700,
        'VSPACING':700,
        'HOVERLAY':0,
        'VOVERLAY':0,
        'CRS':proj_test_etendue,
        'OUTPUT': test_grille_path
    })
 
    test_grille = QgsVectorLayer(test_grille_path, 'test_grille', "ogr")
    project.addMapLayer(test_grille)
    myTestGrille = project.mapLayersByName('test_grille')[0]
    myTestGrille.dataProvider().createSpatialIndex()
 
    # PROCESS CLIP
    _test_clip = myPath + 'clip'
 
    if os.path.isdir(_test_clip) == True:
        shutil.rmtree(_test_clip)
    if os.path.isdir(_test_clip) == False:
        os.mkdir(_test_clip)
 
    test_clip_path = _test_clip + r'\\test_clip.shp'
 
    processing.run("native:clip", {
        'INPUT':test_grille,
        'OVERLAY':test,
        'OUTPUT':test_clip_path
    })
 
    test_clip = QgsVectorLayer(test_clip_path, 'test_clip', "ogr")
    project.addMapLayer(test_clip)
    myTestClip = project.mapLayersByName('test_clip')[0]
    myTestClip.dataProvider().createSpatialIndex()
 
    # PROCESS CENTERLINE
    _test_centerline = myPath + 'centerline'
 
    if os.path.isdir(_test_centerline) == True:
        shutil.rmtree(_test_centerline)
    if os.path.isdir(_test_centerline) == False:
        os.mkdir(_test_centerline)
 
    test_centerline_path = _test_centerline + r'\\test_centerline.shp'
 
    processing.run("Algorithms:Centerlines", {
        'Polygons':test_clip,
        'Method':0,
        'Trim Iterations':10,
        'Simplify':3,
        'Line Spacing':3,
        'Trim Field':'',
        'Simplify Field':'',
        'Densify Field':'',
        'Centerlines':test_centerline_path
    })
 
    test_centerline = QgsVectorLayer(test_centerline_path, 'test_centerline', "ogr")
 
    project.addMapLayer(test_centerline)
    myTestCenterline = project.mapLayersByName('test_centerline')[0]
    myTestCenterline.dataProvider().createSpatialIndex()
 
    # GESTION ORDRE
    root = project.layerTreeRoot()
    root.setHasCustomLayerOrder (True)
 
    order = root.customLayerOrder()
    order.insert(0, order.pop(order.index(myTestCenterline)))
    order.insert(1, order.pop(order.index(myTestClip)))
    order.insert(2, order.pop(order.index(myTest)))
    order.insert(3, order.pop(order.index(myTestGrille)))
    order.insert(4, order.pop(order.index(myTestEtendue)))
    root.setCustomLayerOrder(order)
 
QTimer.singleShot(200, myFunc)