GeoExMachina

Blog géomatique

Tracé Graal (2)

AbracaDijkstra !!!

Suite du Tracé Graal, on continue sur les tracés automatiques de pistes ! Les avancées sont disponibles sur Github, dans le dépôt AutomaTracks.

Adapter un graphe à la création de tracé de route
Malgré mes tentatives de vouloir garder ça simple, ça ne l'est plus beaucoup. En effet le modèle de pistes automatiques que je souhaite réaliser devait respecter un certain nombre de contrainte : pente maximum en long, des coûts en fonction de la pente latérale (surcoût dû au terrassement nécessaire), et des contraintes de direction pour faire en sorte que les pistes soient praticables.

Le paramètre de contrainte de direction est celui qui a entrainé le plus de modification dans la création du graphe. Pour qu'un véhicule puisse progresser sans manoeuvre, il faut que le tracé des pistes respecte un certain rayon de courbure. Un rayon de courbure c'est le rayon du cercle tangeant à un point d'une courbe... vu qu'on travaille sur des polylignes, on a pas de courbe, mais le calcul peut se faire avec les médiatrices de deux lignes successives. Mais pour faire ce calcul il a fallu modifier la nature du graphe.

Auparavant, chaque pixel du MNT était un noeud et ces derniers étaient reliés par des liaisons. Maintenant, un noeud c'est un segment entre deux pixel généré par la méthode de l'article précédent. Pour chaque pixel, on créé 8, 24, 40, 48 noeuds selon le paramètre choisi, et au passage on ne crée pas les noeuds dont la pente excède la pente maximale souhaitée. Chacun de ces noeuds contient l'identifiant du pixel d'où il part et celui du pixel où il va. Ainsi, pour chaque noeud (nd1) on va connaître le pixel d'arrivé et l'on va pouvoir tester la contrainte de direction avec les noeuds (nd2) qui partent de ce même pixel. S'il respecte la contrainte de direction (par défaut j'utilise 20m de rayon de courbure pour un MNT à 5m de résolution), alors on ajoute le noeud (nd2) en liaison au premier noeud (nd1).


Le rayon de courbure varie selon l'azimuth des deux segments


Pour finir, lors de la création de chaque noeud on calcule aussi le coût théorique pour passer du pixel de départ au pixel d'arrivé. La fonction de coût dépend de la pente latérale. A partir de cette pente, on peut calculer l'ampleur du terrassement à faire sur ce segment, c'est le coût considéré dans AutomaTracks.

Déterminer le chemin de moindre coût
Il existe des heuristiques permettant de trouver des solutions de chemin réalisables. Cependant, s'il s'agit de trouver LE chemin de moindre coût entre deux noeuds d'un graphe, il n'existe pas trente-six solutions, il nous faut l'algorithme de Dijkstra.

Celui-ci va partir d'un noeud et il va se propager au fur et à mesure dans le graphe, en retenant le coût le moins élevé pour accéder à chacun des noeuds. Il n'y a pas d'illustration plus parlante que celle-la :



Dans notre cas, chaque noeud a un coût comme attribut, ainsi qu'un coût cumulé qui sert lors de la recherche du chemin. Au début du traitement on initialise le cout cumulé du noeud de départ. Alors l'algorithme calcule le coût d'accès aux noeuds avec lesquels il a une liaison et dont le coût cumulé sera égal au coût cumulé précédent, plus son propre coût. L'algorithme poursuit sa route ainsi jusqu'à tomber sur le noeud d'arrivé.

Utilisation dans AutomaTracks et correction de tracé
Pour nos besoins en Guyane, on a besoin de suivre au maximum les lignes de crêtes afin d'éviter les zones hydromorphes. Par conséquent le réseau de points à relier est issu des polylignes de crêtes. On extrait les pics et les cols remarquables (des points de passages obligés) de ces polylignes, ainsi que les points de départ et d'arrivé de chaque polyligne, afin de garder la cohérence du réseau. Tous nos points ont alors un attribut de ligne L_id et un attribut de point P_id au sein de la ligne. Il est également possible de créer soi-même son réseau de point du moment que les attributs de la couche de points sont respectés.



AutomaTracks va alors réaliser les chemins de moindre coût entre chaque paire de points qui se suivent, en gardant également la cohérence de direction entre deux tracés. Pour cela, il réalise des clips du MNT autour de la zone de calcul. Lors du tracé d'un nouveau chemin de moindre coût, il fera démarrer la propagation du graphe du point de départ, mais aussi depuis les chemins précédemment calculés. Si tel est le cas, AutomaTracks corrige les polylignes concernées pour éviter de laisser des petits morceaux sans cohérence.





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.