GeoExMachina

Blog géomatique

Tracé Graal

Graphe et automatisation

Lors du dernier billet sur Pistecreator (ici), j'ai évoqué l'idée d'un tracé automatique en titre. Ça n'est pas un hasard, c'est un peu ma thématique de travail depuis un moment. Théorie des graphes, algorithme de Dijkstra, chemin de moindre coût, voilà le menu à l'honneur pour tester l'idée d'un tracé automatique. Pour aujourd'hui on va se limiter à un peu de théorie et à convertir un raster en graphe.

Théorie des graphes, kézako ?
La théorie des graphes, c'est une discipline mathématique qui étudie des modèles composés de noeuds et de liaison entre les noeuds. On y reviendra.

Une tâche automatique, ça sous-entend qu'on attend quelque chose de cette tâche et qu'elle le fasse mieux (ou au moins aussi bien et/ou plus rapidement) qu'un être humain. Il faut estimer la qualité du résultat. Bien entendu chacun aura une appréciation différente de la notion de coût, cela peut être un coût économique, un coût environnemental, un mixte des deux... Mais la méthode de calcul va différer selon l'objet de la tâche.
Dans le cas du tracé de piste, l'affaire est compliquée. Il existe des tas de contraintes qu'il faut prendre en compte, comme les pentes (idée fixe de votre serviteur). "Oui certes, les pentes ça joue, mais aujourd'hui on aplanit tout ça au bulldozer"... mais alors ça coûte cher, pour le portefeuille, mais aussi pour l'environnement. Alors autant les prendre en compte ! Prenons notre modèle numérique de terrain (MNT) préféré et traçons nos pistes !


Et ......

Mince, le MNT ne suffit pas, on aura du mal à calculer les coûts d'une piste par pixel, car le pixel en lui-même n'a pas de direction. Donc difficile de calculer des pentes en long, des pentes en travers, d'estimer les coûts de terrassement, le MNT n'est pas approprié tel quel. Bref le MNT est mort, vive le MNT ! Car il va encore nous servir !

Et si on modélisait le MNT en graphes ? Oui, un pixel, un noeud, et des liaisons entre des pixels voisins. Le nombre de pixels voisins est à la discrétion de chacun (et de la précision du MNT aussi). Je vous invite à consulter Jürg Andreas Stückelberger, c'est le plus calé dans le domaine. On peut utiliser les 8 pixels voisins, tout comme les 24, les 40 ou 48 pixels voisins :


Pour chaque pixel, on établit des liaisons avec ses voisins

"Convertir"" un MNT en graphe
Comment on peut faire ça ? Deux petites classes python et une paire de boucle for.
Une classe de graphe :
class Graph():
    #on initialise le graphe vide et on le remplit ensuite
    def __init__(self):
        #une liste pour tous les noeuds du graphe
        self.nodes=[]
    
    #une fonction pour ajouter un noeud avec un simple identifiant
    def add_node(self, id):
        node = Node(id)
        self.nodes.append(node)
Une classe de noeud :
class Node():
    def __init__(self,id):
        #un identifiant
        self.id = id
        #la liste des liaisons
        self.edges= []
    
    #une fonction pour ajouter une liaison entre deux noeuds avec son coût
    def add_edge(self, node, w):
        self.edges.append(node)
        self.weight[node] = w
On change le raster en array :
data = gdal.Open(mnt,1)
proj = data.GetProjection()
scr = data.GetGeoTransform()
resolution = scr[1]
band=data.GetRasterBand(1)
iArray=band.ReadAsArray()
Et pour chaque cellule alors, il suffit de créer un noeud :
G = Graph()
[H,W] = iArray.shape

for i in range(0,H) :
    for j in range(0,W) :
        #on crée un identifiant, qui en plus permet de garder les coordonnées du pixel
        id_node = "x"+str(i)+"y"+str(j)
        #on ajoute la cellule comme noeud
        G.add_node(id_node)
L'identifiant va nous servir à calculer les liaisons possibles. Pour trouver ses voisins, on récupère le x et le y du pixel et on calcule les voisins à partir de vecteur (entre 1 et 48 comme sur le shéma ci-dessus):
#          px  py
shift = [( 0,  0), #0
         (-1,  1), #1
         (-1,  0), #2
         (-1, -1), #3
         ( 0, -1), #4
         ( 1, -1), #5
         ( 1,  0), #6
         ( 1,  1), #7
         ( 0,  1), #8
         (-1,  2), #9
         (-2,  2), #10
         (-2,  1), #11
         (-2,  0), #12
         (-2, -1), #13
         (-2, -2), #14
         (-1, -2), #15
         ( 0, -2), #16
         ( 1, -2), #17
         ( 2, -2), #18
         ( 2, -1), #19
         ( 2,  0), #20
         ( 2,  1), #21
         ( 2,  2), #22
         ( 1,  2), #23
         ( 0,  2), #24
         (-1,  3), #25
         (-2,  3), #26
         (-3,  2), #27
         (-3,  1), #28
         (-3, -1), #29
         (-3, -2), #30
         (-2, -3), #31
         (-1, -3), #32
         ( 1, -3), #33
         ( 2, -3), #34
         ( 3, -2), #35
         ( 3, -1), #36
         ( 3,  1), #37
         ( 3,  2), #38
         ( 2,  3), #39
         ( 1,  3), #40
         (-1,  5), #41
         (-5,  1), #42
         (-5, -1), #43
         (-1, -5), #44
         ( 1, -5), #45
         ( 5, -1), #46
         ( 5,  1), #47
         ( 1,  5)  #48
         ]
On fait donc une deuxième boucle :
for node in G.nodes:
    #supprime le 'x' dans l'identifiant
    id=node.id[1:]
    #sépare les coordonnées x et y en utilisant le 'y' en séparateur
    px,py=id.split('y')
    for sh in shift :
        x_sh,y_sh = sh
        x_end = px + x_sh
        y_end = py + y_sh
        #on vérifie que le noeud peut exister dans l'emprise du raster
        if x_end >= 0 and x_end < H and y_end >= 0 and y_end < W :
            id_end = "x"+str(x_end)+"y"+str(y_end)
            node_end = G.nodes[id_end]
            cout = #là c'est à chacun de mettre sa formule
            #la longueur, à partir des coordonnées;
            #la pente en long à partir des valeur du MNT aux coordonnées des noeuds;
            #une formule de terrassement avec les pentes en travers calculés à 
            #partir des pixels perpendiculaires au vecteur de direction
            node.add_edge(node_end,cout)
Là c'est bon on a un graphe. Nos pixels sont modélisés par des noeuds et on a établi des liaisons entre les noeuds voisins, avec la formule de coût qui va bien. On continuera la prochaine fois pour montrer comment utiliser ce graphe. Pour les plus pressés, le prototype Qgis est sur Github, à poil sans documentation.