GeoExMachina

Blog géomatique

PisteCreator v1.6

Restons dans le droit chemin

Nouvelle option de tracé assisté ! Il s'agît de répondre à une autre situation que celle évoquée dans l'article précédent. Désormais, au lieu de privilégier l'optimisation de la pente en long, on va davantage se préoccuper de la pente en travers.

En effet, pour de la piste de fin de réseau, on ne peut pas mettre les mêmes moyens que pour des pistes structurantes. Pas ou peu de terrassement, on ne compense pas la déclivité du terrain.

Alors gare !... ou le skidder finit comme ça :


Alors on aurait pu chercher la déclivité la plus faible et proposer le segment suivant comme pour la précédente méthode de tracé assisté. Cependant il faut parfois s'écarter cette pente latérale minimale pour se rapprocher d'un bois en particulier. En conséquence, on a opté pour l'affichage d'un cône qui indique là où le tracé respecte les conditions de pentes. La méthode reste la même, on teste les valeurs de pentes sur un arc de cercle de 60° de part et d'autre de l'azimuth du dernier segment validé !

Un exemple reste plus parlant :


Avec le cône on peut englober le plus de points (d'arbres) possible sans risquer l'accident

Par ailleurs si quelqu'un sait pourquoi la fonction setFillColor() ne fait s'applique pas sur les rubberbands, je suis preneur...

Bref voilà ce qui vous attend dans la v1.6 de PisteCreator. En bonus vous remarquerez que l'option de tracé assisté c'est changé en barre volante qu'il est possible de modifier à tout moment. Il y a le mode "cloisonnement" qu'on a vu ici et le mode "échappement" vu auparavant. Le panneau d'option devient également modifiable à tout moment, sans arrêter l'outil d'édition.

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.