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.




PisteCreator v1.5

Vers le tracé automatique ?

Une nouvelle version de PisteCreator disponible et avec elle, un pas vers le tracé automatique. Cette version ajoute une option de tracé assisté ! Encore un changement dans le panneau des options, je ne vais toutefois pas remettre une capture d'écran. Il n'y a qu'une case à cocher et la couleur du tracé assisté en plus !

Attention cependant, ce tracé assisté n'est pas utilisable pour toutes les situations ! En fait elle est plutôt utile pour une situation particulière... Je parle ici des cas où l'on ne peut pas respecter les seuils de pentes voulus, ni en long ni en travers. Dans ces cas là, pas de solutions alternative, il va falloir s'asseoir sur la pente en travers. Pour compenser, il faudra des travaux de terrassements avec remblais et déblais. Je laisse ça aux ingénieurs. L'objectif vers lequel il faut alors tendre, c'est celui de minimiser la distance de route en remblais. Il faut donc monter la pente (ou la descendre) le plus rapidement possible. Cela se traduit par tracer une route qui optimise la pente en long au maximum autorisé. Et c'est exactement ce que fait l'option de tracé assisté !

Et dans ce cas particulier, appuyez sur Entrée en boucle...


Bon... presque, dans les faits il faudra peut-être adapter la distance de calcul du tracé assisté. En effet, pour trouver le meilleur tracé PisteCreator effectue une série de calculs qui dépend :
  • de l'azimuth du dernier segment
  • de la distance maximum tolérée

Pourquoi ? L'azimuth du dernier segment est nécessaire afin de respecter une cohérence dans le tracé de la route. En montée à 10 % de pente, pas question de faire des têtes d'épingle. Faute de calculer directement le rayon de courbure, le calcul de pente va s'effectuer à la distance maximale précisée en option sur un arc de cercle allant jusqu'à 60° maximum de part et d'autre de l'azimuth du dernier segment.


PisteCreator affichera alors le tracé optimal. Cependant, il se peut que ce tracé ne prenne pas en compte un relief intermédiaire entre le point de départ et le point d'arrivé. Il faut alors réduire la distance maximale autorisée. Pour cela, il faut utiliser la touche *. La distance sera alors définie par la longueur qui sépare le dernier point enregistré au curseur de l'utilisateur.


Voilà un exemple d'utilisation ! Après le tracé du premier segment PisteCreator propose le segment suivant en noir. On voit cependant qu'au troisième segment, PisteCreator coupe les isolignes et la pente à mi-chemin du segment est supérieure à 10 %. On adapte alors la distance avec * et on continue à tracer la piste sans problème !

2018 tout en couleurs !

Bonne année 2018 en espérant qu'elle sera fructueuse à tous ! Pour ma part je commence l'année en annonçant qu'il est désormais possible de choisir ses propres couleurs sur PisteCreator. On dira que c'est ma B.A. de début d'année pour les personnes daltoniennes (et pour tous les adeptes de couleurs flashy à souhait).

Le nouveau panneau de configuration ressemble à ça :



Vous retrouverez la release v1.4 sur Github, je referai un upload sur les dépôts Qgis plus tard.