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.
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!