Utilisation de PgRouting dans un calcul de distance🔗

Nous voulons calculer la distance par la route entre deux points. De manière générale, nous voulons calculer plus d'une distance: il faudra le faire pour un tableau de données qui contient près de 10000 lignes. Ce tableau ne comporte pas de coordonnées géographiques mais uniquement des codes INSEE.

C'est un cas assez général: on souhaite faire un calcul de distance par la route entre deux communes.

Pour réaliser ce calcul, nous allons utiliser les outils suivants:

• PostgreSQL pour stocker les données attributaires
• PostGIS pour stocker les données géographiques comme le réseau routier et les points de départ et d'arrivée
• pgRouting qui est l'extension de calcul

pgRouting implémente plusieurs méthodes pour trouver le chemin le plus court entre deux points d'un réseau routier.

Ce document a pour objectif de présenter une méthode complète et détaillée du comment-faire . Le but est d'être réutilisable dans d'autres conditions.

• Postgresql v8.4 (administrateur = compte postgres)
• PostGIS 1.5
• pgrouting 1.05
• Une base de données spatiale déjà installée (nommée GEOBASE_LOCALE dans le document)
• la BDROUTE500

Vous devez déjà disposer d'une DB PostGIS fonctionnelle ! Pour l'instant, pgrouting n'est pas intégré comme un paquet Debian. Vous devez-donc le télécharger depuis son site web et le compiler pour enfin l'installer. Ensuite, il faudra importer les fonctions dans votre base de données spatiale.

Les commandes qui suivent s'appliquent au système d'exploitation Debian en version testing (nom de la distribution Wheezy). Ce document essaye de se focaliser sur la partie données. Si vous voulez consultez les instructions d'installation de pgRouting sous Debian, vous pouvez lire ce qui suit…

Pour compiler pgRouting, vous aurez besoin de cMake et d'un patch:

# aptitude install cmake libcgal-dev postgresql-server-dev-8.4 libboost-graph-dev 
$ wget http://download.osgeo.org/pgrouting/source/pgrouting-1.05.tar.gz
$ tar xzf ./pgrouting-1.05.tar.gz
$ cd pgrouting-1.05
$ cmake -DWITH_DD=ON .
$ patch -p1 -i ~/pgrouting_debian_wheezy_boost1.48.patch
$ make
# make install

Normalement, l'installation a déposé des scripts SQL qui permettent de charger les fonctions pgrouting dans une base de données PostGIS. Nous chargeons également les fonctions liées aux fonctionnalités Drive-Distance.

$ su postgres
$ psql -U postgres -f /usr/share/postlbs/routing_core.sql GEOBASE_LOCALE
$ psql -U postgres -f /usr/share/postlbs/routing_core_wrappers.sql GEOBASE_LOCALE
$ psql -U postgres -f /usr/share/postlbs/routing_topology.sql GEOBASE_LOCALE
$ psql -U postgres -f /usr/share/postlbs/routing_dd.sql GEOBASE_LOCALE
$ psql -U postgres -f /usr/share/postlbs/routing_dd_wrappers.sql GEOBASE_LOCALE
$ exit

Nous allons charger plusieurs couches de la BDROUTE500 dans PostGIS. Vu le poids des fichiers, on peut charger l'ensemble des données. Nous allons utiliser ogr2ogr pour intégrer ces couches au format MapInfo.

Deux couches nous intéressent:

• Celle du réseau (c'est la base)
• Celle des noeuds routiers en sachant que nous désirons uniquement récupérer les noeuds correspondant aux centre des communes.

Voici les séquences ogr2ogr pour intégrer ces données dans notre base de données spatiales. Vous noterez qu'à chaque fois, nous sommes obligés de forcer la projection (options -s_srs et -a_srs) car la gestion des couches MapInfo met en oeuvre d'autres références de référentiels. Si on laisse les options par défaut, il n'y a pas de correspondances avec les codes EPSG (même si les données sont dans le bon référentiel). Cela a des conséquences sur les requêtes spatiales PostGIS qui imposent de travailler dans le même référentiel pour toutes les couches.

Intégration des routes:

ogr2ogr -f "PostgreSQL" PG:"host=localhost user=sid_admin dbname=GEOBASE_LOCALE password=mypassword" ~/GEOBASE/GB_REF/BDROUTE500/RESEAU_ROUTIER/N_ROUTE_RT5_000.TAB -nln route -a_srs EPSG:2154 -s_srs EPSG:2154

Intégration des points de départ/d'arrivée des communes:

ogr2ogr -s_srs EPSG:2154 -a_srs EPSG:2154 -f "PostgreSQL" PG:"host=localhost user=sid_admin dbname=GEOBASE_LOCALE password=ddaf49" ~/PT_COMMUNE.TAB -nln nd_com_049

Pour les besoins de croisement, il nous faut la couche des emprises de commune. Nous allons utiliser la BDCARTO. Il nous faut également l'importer dans GEOBASE_LOCALE:

ogr2ogr -f "PostgreSQL" PG:"host=localhost user=sid_admin dbname=GEOBASE_LOCALE password=mypassword" ~/GEOBASE/GB_REF/BDCARTO/ADMINISTRATIF/N_COMMUNE_BDC_049.TAB -nln communes -a_srs EPSG:2154 -s_srs EPSG:2154

Avant de nous lancer sur la table complète des routes de France, il est bon de fabriquer un extrait de cette couche pour réduire les temps de calcul.

CREATE TABLE route_049 As SELECT route.* FROM route, communes WHERE ST_Intersects(route.wkb_geometry, communes.wkb_geometry);
SELECT Populate_Geometry_Columns();

Avant de lancer les calculs avec pgRouting, nous devons transformer la couche des routes en réseau topologique car pgrouting a besoin d'informations supplémentaires, notamment les points d'arrivées et de départ. De plus, nous aurons besoin de nommer correctement ces points en fonction d'un code INSEE.

Pour construire un réseau topologique, il faut d'abord créer des colonnes supplémentaires:

ALTER TABLE route_049 ADD COLUMN "source" integer;
ALTER TABLE route_049 ADD COLUMN "target" integer;

Ensuite, nous allons utiliser une fonction de pgrouting pour créer une couche supplémentaire de géométrie qui sera utilisée lors des calculs. C'est la fonction assign_vertex_id:

SELECT assign_vertex_id('route_049', 0.00001, 'wkb_geometry', 'id_route50');

Au final, nous obtenons des numéros d'identification de points dans les colonnes source et target. Ensuite, nous avons également une table supplémentaire nommée vertices_tmp qui contient des points avec des identifiants.

Le nommage en fonction du code INSEE est une simple requête géographique avec la couche des communes que nous avons intégré au départ:

SELECT communes.insee_com::integer, nd_com_049.wkb_geometry FROM nd_com_049 INNER JOIN communes ON ST_Intersects(nd_com_049.wkb_geometry,communes.wkb_geometry)

Pour conserver cette table, réalisons une table de liaison de la forme avec l'instruction suivante (la structure est simplissime):

CREATE TABLE vertices_insee As
  (SELECT a.id, b.insee FROM
   vertices_tmp As a
   JOIN 
   (SELECT communes.insee_com As insee, nd_com_049.wkb_geometry
     FROM nd_com_049 
     INNER JOIN communes ON ST_Intersects(nd_com_049.wkb_geometry,communes.wkb_geometry)
   ) As b
   ON a.the_geom=b.wkb_geometry
  );

Ainsi, lorsque nous voudrons lancer nos calculs finaux, il nous faudra juste utiliser cette table intermédiaire pour faire le lien entre un code INSEE et un identifiant de point de la table vertices_tmp. On aurait également pu ajouter une colonne code_INSEE dans la table vertices_tmp ou faire des sous-requêtes dans les requêtes à suivre. Néanmoins, cette table nous permettra de faciliter les explications qui suivent.

Maintenant que tout est prêt, faisons un petit essai: essayons de calculer la distance entre deux points:

SELECT * FROM shortest_path('
     SELECT id_route50 as id,
	      source::integer,
	      target::integer,
	      longueur_t::double precision as cost
	     FROM route_049',
     753, 1926, false, false);

Le résultat est présent sous forme de tableau:

SELECT sum(b.cost) As longueur FROM (
SELECT * FROM shortest_path('
	   SELECT id_route50 as id,
		    source::integer,
		    target::integer,
		    longueur_t::double precision as cost
		   FROM route_049',
	   753, 1930, false, false)
) As b;

Nous voulons calculer les distances kilométriques avec pgrouting entre toutes les communes de notre département. Le département d'étude est le Maine-et-Loire, il comprend 363 communes (donc 363 codes INSEE). Cela nous donne 363^2 possibilités (soit 131769 combinaisons possibles).

Nous devons déjà réaliser une table comprenant toutes ces lignes:

CREATE TABLE distance_departement
(
  id serial,
  depart character(5),
  id_depart integer,
  arrive character(5),
  id_arrive integer,
  distance real
);

Ensuite, nous allons la remplir avec les valeurs qui vont bien (les codes INSEE sont directement issus de la table communes). Une simple requête mal formée suffit pour combler notre besoin. Il faut ensuite également ajouter les données des identifiants des points pour faire le lien entre les codes INSEE et les identifiants de points de la table vertices_tmp. Nous allons, pour ce faire, chercher les données de la table spécialement créée pour l'occasion: vertices_insee.

-- Création des lignes de code INSEE
INSERT INTO distance_departement (depart,arrive) SELECT a.insee_com,b.insee_com FROM communes As a, communes As b;
-- Lien avec les identifiants de la table vertices_tmp en départ
UPDATE distance_departement SET id_depart = b.id 
FROM vertices_insee As b WHERE b.insee = distance_departement.depart
-- Lien avec les identifiants de la table vertices_tmp en départ
UPDATE distance_departement SET id_arrive = b.id 
FROM vertices_insee As b WHERE b.insee = distance_departement.arrive

Nous devons maintenant calculer la distance pour chaque ligne de la table distance_departement. Nous allons donc réaliser une petite requête UPDATE assez classique:

UPDATE distance_departement
SET distance = (SELECT sum(b.cost) As longueur FROM (
 SELECT * FROM shortest_path('
	    SELECT id_route50 as id,
		     source::integer,
		     target::integer,
		     longueur_t::double precision as cost
		    FROM route_049',
	    id_depart, id_arrive, false, false)
 ) As b)
WHERE id_depart IS NOT NULL AND id_arrive IS NOT NULL AND id_depart != id_arrive;

Attention, cette requête prend un temps non négligeable… Comptez environ 2-3 heures pour un département ! De plus, elle consomme un sacré paquet de mémoire. L'idéal est donc de la découper:

UPDATE distance_departement
SET distance = (SELECT sum(b.cost) As longueur FROM (
 SELECT * FROM shortest_path('
	    SELECT id_route50 as id,
		     source::integer,
		     target::integer,
		     longueur_t::double precision as cost
		    FROM route_049',
	    id_depart, id_arrive, false, false)
 ) As b)
WHERE depart IN (SELECT insee FROM vertices_insee 
WHERE insee::integer > 49050 and insee::integer < 49100
ORDER BY INSEE) AND id_arrive IS NOT NULL AND id_depart != id_arrive

Vous devez utiliser une version récente de Calc car vous aurez plus de 65535 lignes à gérer. Dans ce cas, je vous préconise d'utiliser celle de LibreOffice.

Grâce à ce document et à un peu de logiciels libres, nous voilà en mesure d'effectuer des calculs de distance entre deux communes ou deux points sur un réseau routier. C'est un problème courant dont la résolution dépend en partie de la qualité de l'information qui est utilisée.

Dans notre cas, nous avons utilisé la BDROUTE 500 ainsi qu'une approximation des points qui constituent les centres des communes. On pourrait imaginer un réseau plus complet et surtout, une base de données plus précise des points de départ et d'arrivée sur lesquels on souhaite faire des calculs.

Si votre structure est une administration, vous pouvez, à l'aide de ces données et de la connaissance des flux de citoyens qui viennent vous consulter directement, faire une estimation un peu plus précise du bilan carbone engendré par cette partie spécifique.

Si vous réalisez des études sur le domaine des déplacements, vous pouvez chiffrer plus exactement quelles sont les distances parcourues et ce, pour un volume conséquent d'informations.

Un tas d'autres utilisations et couplages sont possibles. Maintenant que vous avez les outils et la méthode, vous voilà bien servis pour répondre à vos nombreuses interrogations sur les calculs de distance.

• Beginners guide to PgRouting
• Workshop PgRouting
• Pgrouting sur le site de l'université Stanford
• Pgrouting chez cartoweb