lundi 14 mars 2016

ArcGIS: Créer un index des voies pour un atlas communal

On travaille actuellement à créer une base de données du réseau routier de Tahiti. Les communes de Tahiti sont des partenaires importants, avec lesquels on valide les tracés des voies et leurs noms. En retour, nous fournissons des shapefiles du réseau et un atlas en pdf de la voirie pour chaque commune.

Exemple de carte de l'atlas d'une commune


Pour fixer les choses, la base de données de travail du réseau routier est au format géodatabase fichier ESRI. On travaille avec ArcGIS 10.x et le module NetworkAnalyst pour les analyses les plus poussés. Les atlas des voies sont réalisés avec les outils des "pages dynamiques" d'ArcGIS. Pour que l'atlas présente bien, on réalise aussi une page de couverture. On a voulu aussi ajouter un index des voies: une liste alphabétique des voies avec les carreaux dans lesquels elles figurent. Ce dernier point a présenté une certaine difficulté. ArcGIS ne propose pas d'outil permettant de créer directement une telle table. On a pensé pouvoir le réaliser sous Excel à partir d'un export des voies croisées avec les carreaux. Mais ce n'est pas si simple qu'on le pensait!
La colonne index donne la liste des cartes de chaque voie 

Pour faciliter la production de l'index, j'ai réalisé un script python pour construire cette table. Bien sûr, le script est ensuite intégré dans une boîte à outils dédiée, pour plus de facilité. L'idée est de partir de la classe d'entités des segments des voies (ou des itinéraires), de la classe d'entités des index des cartes (les carreaux numérotés des emprise des cartes) et de construire une table à deux attributs:
  • l'identifiant des voies
  • la liste concaténée des cartes où figure la voie
On fera ensuite des jointures sur d'autres tables puis on retravaillera le résultat dans Excel pour un rendu de qualité. Au passage, on a essayé d'utiliser l'outil de rapport fourni avec ArcGIS, mais il est beaucoup moins souple d'Excel.

L'algorithme du script Python repose sur une boucle de parcours de toutes les voies. toutes les voies qui ont un nom ou un numéro doivent figurer dans l'index. Les voies sans nom n'y figurent pas. On sélectionne tous les segments d'une voie et on fait une intersection par emplacement avec la couche des index de cartes (les carreaux numérotés des emprises). On concatène la liste des numéros des cartes dans un seul attribut (appelé "index"...).

L'attribut itinéraire correspond au identifiant des voies.
Les valeurs mini et maxi permettent de limiter les rues à faire figurer dans l'index (ceux d'une même commune).
L'attribut de carroyage est celui des noms (ou numéros) des cartes.
Les deux derniers paramètres sont là pour créer la table d'index.
Le script Python de cet outils peut être téléchargé ici: AtlasIndex.py



dimanche 21 février 2016

PostGIS: orienter un symbol ponctuel le long d'une route

On a levé les emplacements des portails sur le réseau routier en même temps qu'on a constitué notre base de données des routes. Ces portails représentent des entrées de lotissements sécurisés, des accès à des zones réglementées, etc. Ils empêchent la circulation des véhicules et sont donc importants à représenter sur une carte.

Les portails sont dans une couche de points et les routes dans une couche de lignes. L'image suivante montre les deux couches:
On se pose un problème de représentation des portails: on voudrait un symbole qui soit toujours perpendiculaire à la voie. Pour cela, il faudrait connaître l'orientation de la route au niveau de chaque point représentant un portail... Les deux couches étant indépendantes, ce n'est pas très simple!

Les deux couches sont chargées dans une base de données PostgreSQL/PostGIS. La couche de point s'appelle 'portail' et la couche des routes 'segment'. Avec la requête suivante, on trouve le segment (s.*) correspondant à chaque chaque portail (p.*):

select p.*, s.*
from portail p
join segment s on s.id = (select id
        from segment
        order by ST_Distance(geom, p.geom)
        limit 1) ;


Cependant, la couche 'segment' est topologique: chaque segment va d'un carrefour à un autre (pour faire simple). Un segment peut avoir beaucoup de sommets (vertex) et une orientation très variable. Il faut donc "descendre" plus finement dans le segment pour regarder l'orientation au niveau du portail.

Deux fonctions PostGIS vont nous aider:
  • La fonction ST_Line_Locate_Point(s.geom, p.geom) donne la distance (en pourcentage) du portail (p.geom) par rapport au début du segment (s.geom).  
  • La fonction ST_Line_Interpolate_Point(s.geom, distance) fait le travail inverse: elle fournit un point sur le segment 's.geom' à partir de la distance (en pourcentage) du début.
L'idée consiste à s'inspirer du calcul de la pente en un point: on cherche les coordonnées d'un point sur le segment un peu avant (5 m) le portail et d'un point un peu après (toujours 5m). On calculera l'azimuth entre ces deux points pour l'assimiler à la direction au niveau du portail.

La distance de 5 m est retenue en considérant l'espacement des sommets de notre réseau routier, car il est rarement inférieur à cette distance. Il y a quand même des cas où le portail est à moins de 5m du début ou de la fin du segment... Il ne faudra pas oublier de traiter ces cas! 

Il y a une petite conversion à effectuer pour transformer la distance de 5m en relatif à la longueur du segment, puisque la fonction ST_Line_Locate_Point et son corolaire utilise une distance en pourcentage et non la longueur réelle du segment. Notre réseau routier est en coordonnées géographique (WGS84), il faut aussi convertir la distance du segment en mètre... Le calcul suivant nous sort d'embarras:  5/ST_Length(ST_Transform(s.geom, 3395)) 
ST_Transform(s.geom, 3395) projette en "World Mercator" (EPSG:3395) le segment 's.geom'.

Enfin la fonction ST_Azimuth permet de calculer l'angle par rapport au nord entre les deux points entourant le portail et on ajoute 90° pour avoir la perpendiculaire à la route.

J'ai mis tout ça dans une fonction 'point_direction' téléchargeable: fonction_point_direction.sql

Je peux maintenant calculer la colonne 'angle' de la couche 'portail' avec la requête de mise à jour suivante:

update portail p
set angle = (select point_direction(p.geom, s.geom, 'true')
        from segment s
        where s.id = (select id
                from segment
                order by ST_Distance( geom, p.geom)
                limit 1)
        ) ;


Enfin, j'arrive à représenter les portails avec un symbole orienté correctement par rapport à la route!

  








lundi 15 février 2016

PostGIS et pgrouting: Comment ordonner les segments d'une rue ?

On a chargé dans PostGIS une table des segments des rues (table "segment"). Cette table a été construite sous ArcGIS et transférer dans PostGIS avec FME. Le réseau que l'on a créé est orienté: le sens positif correspond au sens de circulation. Lorsque les voies sont à double sens, une convention est utilisée pour orienter les segments. Ainsi, tous les segments d'une même voie (rue ou route) sont dans le même sens.

Dans PostGIS, on utilise l'extension pgRouting pour travailler avec ce réseau. La première étape est de créer la topologie avec la fonction "pgr_CreateTopology". Cette fonction ajoute deux colonnes "source" et "target" à la table "segment" (celle qui contient la géométrie des segments) et elle crée un nouvelle table "segment_vertices_pgr" (le préfixe "segment" provient du nom de notre table). C'est une table de points: un champ géométrie de type point et un identifiant, les autres champs ne sont pas importants. Ces points sont les jonctions entre les segments, c'est-à-dire les nœuds du réseau topologique. La fonction "pgr_CreateTopology" ajoute les identifiants des nœuds de début et de fin de chaque segment dans les champs "source" et "target".

Dans la table des segments, j'ai aussi des informations sur les codes des rues. Une rue est généralement composée de plusieurs segments. Par exemple, je fais une requête pour retrouver tous les segments de la rue 380027 (c'est son code: son nom n'a pas d'importance ici):

select id, source, target
from route.segment
where code_rue = 380027
order by id ;


Ce qui me renvoie le tableau suivant:


id source target
23348 14300 4332
23359 14331 14299
23419 14301 14362
23420 14299 14368
23421 14362 14300
23489 14368 14301

Avec QGIS, voila ce que ça donne:


En blanc, on a représenté la table "segment" avec l'identifiant comme label. On a aussi ajouté la table "segment_vertices_pgr" dont les identifiant sont soulignés. La rue 380027 est en jaune sur l'image et commence à gauche. Le premier segment est le 23359, puis 23420, 23489, 23419, 23421 pour finir par le segment 23348.

On voit que les segments ramenés par la requête ne sont pas dans l'ordre d'affichage. Sur l'image, les numéros des segments (les "id" de la requête) ne se suivent pas. Il n'y a pas d'ordre non plus dans les identifiants des nœuds. Le "order by" de la requête ne sert à rien!

On sait que les champs "source" et "target" donnent une information sur le sens du segment: le nœud "target" du premier segment de la rue est le nœud "source" du second segment, etc. Les segments sont chaînés grâce à ces deux champs, mais quelle requête pourrait les lister dans l'ordre d'affichage ?

Ce problème de chaînage est classique en informatique. On utilise souvent un algorithme récursif de parcourir d'arbre. La documentation de PostgreSQL fourni un mécanisme de pré-requête avec "WITH" et une variante intéressante pour notre cas "WITH RECURSIVE". Cette pré-requête a pour intérêt d'éviter une sous-requête (syntaxe lourde et peu performant). Le fonctionnement récursif n'est pas très  bien expliqué dans la documentation de PostgreSQL. Par contre, le blog de Tran Xuan Truong donne un exemple d'utilisation de "WITH RECURSIVE" pour parcourir un arbre.

Pour qu'un algorithme récursif "fonctionne", il y a trois conditions à respecter: une condition de départ, une condition de propagation et une condition d'arrêt. Dans notre cas, la condition de départ est donnée par le premier segment: c'est celui dont le nœud "source" ne se retrouve pas dans les "target" des segments de la rue (il n'est chaîné à aucun segment de la rue):

select s.id, s.source, s.target
from route.segment s
where s.code_rue = 380027
and s.source not in (select s.target
    from route.segment s
    where code_rue = 380027) ;

 
Ce qui renvoie:

id source target
23359 14331 14299

La condition de propagation est donnée par le chaînage: le noeud "source" du segment suivant est le noeud "target" du segment précédent. La condition d'arrêt est ici assez simple puisque la liste des segments de la rue est finie (5 segments).

Ces trois conditions sont mis en forme dans la requête SQL de la manière suivante:

WITH RECURSIVE ordered_list(id, source, target, seq) as (
    select id, source, target, 1
    from route.segment
    where code_rue = 380027
    and source not in (select target
        from route.segment
        where code_rue = 380027)
    UNION ALL
    select s.id, s.source, s.target, ol.seq + 1
    from route.segment s
    join ordered_list ol on s.source = ol.target
    where code_rue = 380027
    )
select *
from ordered_list ;



L'ordre "WITH RECURSIVE" crée une table temporaire, que j'appelle "ordered_list". Dans la première partie de la définition, on retrouve la condition de départ. Ensuite, "UNION ALL" permet d'ajouter la condition de propagation: un "select" sur les tables "segment" et "ordered_list" (la récursivité s'applique ici) et une jointure entre les nœuds "target" et "source". On peut remarquer que j'ai ajouté un attribut "seq" pour suivre le parcours.

Le résultat est le suivant:


idsourcetargetseq
2335914331142991
2342014299143682
2348914368143013
2341914301143624
2342114362143005
233481430043326

Grâce à cette ordre "WITH RECURSIVE" on peut parcourir tout un réseau dans l'ordre. Cette extension du langage SQL est particulièrement intéressante pour les réseaux où ce type de problème se rencontre fréquemment.
 




vendredi 5 février 2016

Définir la tolérance de la topologie dans ArcGIS

La topologie d'ArcGIS est un super outil, bien utile pour nettoyer les données et les structurer géométriquement. Il y a deux paramètres spécifiques à la topologie: la précision et la tolérance (ou tolérance d'agrégat, traduction de "cluster tolerance").

La précision est facile à comprendre: c'est une grille sur laquelle sont calées (donc modifiées) les coordonnées en entrée. Elle ne doit pas être trop grande pour affecter la précision réelle des données, mais ce n'est pas la peine non plus de fixer une précision trop importante. Le conseil généralement donné est de fixer la précision au 1/10 de la tolérance. Ok, mais la tolérance, c'est quoi ?

La documentation d'ArcGIS est un peu vague sur la tolérance d'agrégat. Comme pour la précision, il s'agit d'une grille, donc la tolérance est une largeur en x et en y. Le problème est que, pour tout cartographe, la tolérance n'est pas comprise comme ça. Elle dépend de l'échelle nominale de la couche, c'est à dire l'échelle "normale" pour laquelle les données sont produites. C'est une notion importante en SIG qu'il faut toujours avoir à l'esprit.

Je vais prendre l'exemple d'une couche de parcelles du cadastre établie au 1/1.000. On considère généralement que la tolérance sur les données à cette échelle est de 2/10 mm papier, soit 20 cm terrain. Comment fixer la tolérance d'agrégat ? Les différents tests que j'ai mené me conduise à utiliser la formule suivante: tolérance d'agrégat = tolérance des données / (2 * racine(2) ). Pour la couche parcelle à l'échelle 1/1.000, la tolérance des données est de 20 cm, donc la tolérance d'agrégat est de 7 cm. Le terme racine(2) vient de ce que la tolérance d'agrégat est une grille et notre tolérance est une distance entre deux points: il faut appliquer Pythagore pour faire apparaître ce facteur.
La tolérance d'agrégat en x et y. La tolérance des données en diagonale.
 Dans l'exemple suivant, j'ai créé une couche de polygones et j'ai fixé la tolérance d'agrégat à 7 cm. Le cercle rouge matérialise la tolérance des données (rayon de 20 cm).
Les deux rectangles se touchent au centre du cercle. La pointe du triangle est en dehors du cercle de précision. La validation topologique ne change rien, ce qui est normal.

Attention: la topologie modifie aussi tous les vertex des formes. Dans les exemples, j'ai espacé les vertex des formes de plus de 20 cm. Sinon, ils seraient modifiés par la validation topologique.

Maintenant, je déplace légèrement la pointe du triangle pour qu'elle entre dans le cercle rouge:
Et je fais une validation topologique: Les trois sommets sont déplacés pour ne faire qu'il seul point.


Un autre cas de figure:
Cette fois, c'est le rectangle de droite que je déplace. Le sommet de celui de gauche est au centre du cercle et la pointe du triangle est à l'extérieur du cercle.

La validation topologique rassemble les sommets des deux rectangles, sans modifier le triangle.

Un dernier exemple:
Le sommet du rectangle de gauche est au centre du cercle rouge, les 2 autres figures sont à l'extérieur.
La validation topologique ne change rien (l'image ne change pas, je ne la répète pas!).

En conclusion, on a une règle pour fixer la tolérance d'agrégat de la topologie ArcGIS en fonction de l'échelle nominale des données:
  • Echelle: 1/1.000, tolérance d'agrégat: 0,07m
  • Echelle: 1/2.000, tolérance d'agrégat: 0,14m
  • Echelle: 1/5.000, tolérance d'agrégat: 0,35m
  • Etc.
Il faut être conscient que tous les points à une distance entre eux inférieure à la tolérance seront modifiés: les sommets comme les vertex. Cela peut dérouter au début, mais il faut être conscient que cela n'aura pas d'impact visuel sur les données, tant qu'elles sont utilisées à des échelles réalistes! Ceux qui zoome trop (si, si, ça se fait!) ne verront plus les petits trous ou superpositions, qui donne l'impression que la données est pleine d'erreur.
Autre avantage: les données sont allégées. Plus de points utile! La topologie fait un sacré ménage!


jeudi 12 novembre 2015

Semaine de l'innovation publique 2015

Mon équipe et moi-même avons eu le privilège d'être sélectionné parmi les applications "innovantes" de l'administration.
Cette manifestation s'inscrivait dans le cadre de la semaine de l'innovation publique, relayée en Polynésie par le Haut-Commissariat et la Présidence. La semaine se résumait à une journée, le 15, où des présentations permettaient au services de se retrouver mais aussi le grand public, puisque la présidence était ouverte à cette occasion.
Pour élargir l'audience de la présentation, je l'ai chargé sur SlideShare.net et vous pouvez la consulter en version encapsulée ici:



mardi 13 octobre 2015

Te Fenu@ dans un site

On nous demande souvent comment intégrer Te Fenu@ dans un site web, à la manière de google map. Pour cela, il y a une URL dédiée : https://www.tefenua.gov.pf/tefenua/embed.html

Cette URL peut être intégrée dans une page HTML avec une balise <iframe> comme ci-dessous :

<iframe src="https://www.tefenua.gov.pf/tefenua/embed.html?zoom=8" width="600" height="300"></iframe>

Ce qui donne ceci:



Pour l'instant cette URL accepte les paramètres suivants :
- zoom: un entier de 0 à 18. Plus le nombre est grand, plus l'échelle est grande. zoom=18 correspond à l'échelle max, soit 1:1.000.
- center: le centre de la carte en longitude et latitude.
- point[]: un tableau de point avec un marqueur sur chaque point.

L'exemple suivant centre la carte sur le point (long=-149.6107, lat=-17.6088) à l'échelle 1:57.000 (zoom=12) et affiche 2 points repérés par un marqueur (point[0] et point[1]) :

https://www.tefenua.gov.pf/tefenua/embed.html?center=-149.6107,-17.6088&zoom=12&point[0]=-149.6107,-17.6088&point[1]=-149.6087,-17.62

Pour le résultat suivant:

 

Dernière petite idée que je voudrais vous glisser à l'oreille: il est aussi possible de faire un lien dans un mail, par exemple sur une position dans Te Fenu@ :

https://www.tefenua.gov.pf/tefenua/?center=-149.6107,-17.6088&zoom=12&point[0]=-149.6107,-17.6088&point[1]=-149.6087,-17.62

L'exemple est le même que précédemment, mais cette fois, il n'y pas "embed.html" dans l'url. Les autres paramètres restent les mêmes.