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!