Annexe A. PgRouting pour le calcul d'itinéraire

Table des matières

A.1. Installation sous Windows
A.2. Installation sous Ubuntu Dapper/Edgy
A.2.1. GAUL
A.2.2. BGL
A.2.3. CGAL
A.2.4. PgRouting
A.3. Installation sous Mac OS X
A.3.1. FINK
A.3.2. GAUL
A.3.3. BGL
A.3.4. CGAL
A.3.5. PgRouting
A.4. Chargement des fonctionnalités de PgRouting
A.5. Importation d'un shapefile concernant des tronçons
A.6. Obtention des noeuds du réseau
A.7. Fonctionnalité shortest_path_astar()
A.7.1. Exemple pour les noeuds 38 et 48.
A.7.2. Démo en ligne avec MapServer
A.7.3. Tester soi-même la démo avec MapServer.
A.7.4. Tester sur un jeu de données réelles: jeu de tests GEOROUTE IGN
A.8. Fonctionnalité TSP()
A.8.1. Exemple
A.8.2. Programme en C pour les appels successifs à shortest_path_astar()
A.8.3. Limites du programme
A.9. Fonctionnalités shortest_path_astar_as_geometry_internal_id_directed() et tsp_astar_as_geometry_internal_id_directed()
A.9.1. Importation d'un jeu de données NavTeq
A.9.2. Noeuds du réseau et direction pour le routage
A.9.3. Modifications nécessaires sur la table streets_edges
A.9.4. Exemple avec shortest_path_astar_as_geometry_internal_id_directed()
A.9.5. Exemple avec tsp_astar_as_geometry_internal_id_directed()

Cité en exemple dans l'Avant-propos de la documentation, PgRouting est la librairie développée en étroite collaboration par Camptocamp (France/Suisse) et Orkney (Japon). PgRouting fait partie du projet PostLBS. Le site de PgRouting est http://www.postlbs.org/postlbs-cms/en/project.Gérald FENOY a proposé une traduction de la documentation de pgrouting disponible à http://www.postgis.fr/book/print/360

Je propose ici de montrer un exemple d'utilisation de la fonctionnalité shortest_path_astar() et la fonctionnalité tsp(). Comme l'exemple ici proposé est très minimaliste, je ne prends pas en considération l'usage des index! Merci d'en prendre note.

A.1. Installation sous Windows

Pour PostgreSQL 8.2.X, il faut se procurer le fichier zippé http://www.postlbs.org/postlbs-cms/files/downloads/pgRouting-0.9.8-0_win32.zip. Une fois dézippé, il suffira de le copier-coller dans votre arborescence des répertoires (share, lib, doc) de PostgreSQL.

A.2. Installation sous Ubuntu Dapper/Edgy

Conformément à la documentation sous Ubuntu, l'installation n'est pas des plus compliqués.Rapidement voyons les diverses étapes de l'installation

A.2.1. GAUL

Il faut se rendre à http://prdownloads.sourceforge.net/gaul/gaul-devel-0.1849-0.tar.gz?download pour télécharger gaul-devel-0.1849-0.tar.bz2. Ensuite

tar xzf gaul-devel-0.1849-0.tar.gz
cd gaul-devel-0.1849-0
./configure --enable-slang=no
make
make install

A.2.2. BGL

La commande suivante suffira amplement

apt-get install libboost-graph-dev libboost-graph1.33.1

A.2.3. CGAL

Il faut télécharger et l'installer

wget ftp://ftp.mpi-sb.mpg.de/pub/outgoing/CGAL/CGAL-3.2.1.tar.gz
tar xvzf CGAL-3.2.1.tar.gz 
cd CGAL-3.2.1
./install_cgal --prefix=/usr/local/cgal --with-boost=n --without-autofind -ni /usr/bin/g++

Puisque la librairie routing.so que nous allons installée sera lieé par la suite à libCGAL.so, prenons déjà les devants. En effet, lors de l'appel de la librairie routing.so par le biais des fichiers SQL routing.sql et routing_postgis.sql, cette librairie a besoin de savoir où se trouve libCGAL.so

echo /usr/local/cgal/lib/i686_Linux-2.6_g++-4.0.3/ >> /etc/ld.so.conf
ldconfig

A.2.4. PgRouting

Nous ferons tous simplement

wget http://www.postlbs.org/postlbs-cms/files/downloads/pgRouting-0.9.9.tgz
tar xvzf pgRouting-0.9.9.tgz
cd routing/
export CGAL_MAKEFILE=/usr/local/cgal/make/makefile_i686_Linux-2.6_g++-4.0.3
./configure --with-cgal=/usr/local/cgal --with-gaul=/usr/local/
make
make install

A.3. Installation sous Mac OS X

L'installation sous cet OS nécessite quelques subtilités, fruit de test du travail de Nicolas RIBOT, que je salue ici au passage ;).

A.3.1. FINK

Il vous faut installer Fink. Puis depuis Fink, faites "Fink -> preferences -> Fink -> Check: Use unstable packages"

A.3.2. GAUL

Télécharger la dernière version à http://prdownloads.sourceforge.net/gaul/gaul-devel-0.1849-0.tar.gz?download . Puis faites

export CFLAGS=
./configure --enable-ccoptim=no  --enable-slang=no
make
sudo make install

A.3.3. BGL

Télécharger BOOST libs: www.boost.org, version 1.33.1. Excéutez ensuite les commandes suivantes

./configure
make

Pensez à aller prendre un café! Puis faites

make install

Un deuxième café ne serait pas de trop.Il faut maintenant compiler bjam. Pour celà rendez-vous dans tools/build/jam_src. Faites ensuite

./build.sh

Il faut créer un lien symbolique depuis tools/build/jam_src/bin.macosxx86/bjam dans un folder contenu dans le path

cd /usr/bin; ln -s <path_to_bjam>

Puis exécutez

bjam install

Celà copiera les headers dans /usr/share/boost_1_33_1/. Exécutez ensuite

bjam stage 

qui copiera les les lib dans un seul endroit.

A.3.4. CGAL

Télécharger le DMG pour Mac http://www.cgal.org/cgi-bin/cgal_download.pl Installez-le en faisant

CXX=/usr/bin/g++ ./install_cgal -ni --BOOST_INCL_DIR /usr/local/include/boost-1_33_1 \
--without-autofind --BOOST_LIB_DIR /Users/nicolas/download/boost_1_33_1/stage -with-BOOST

A.3.5. PgRouting

Mettre a jour les chemins vers CGAL, BOOST, GAUL. créer un lien symbolique pour BOOST dans /usr/local:

ln -s boost-1_33_1/boost .

Créer un lien symbolique dans /usr/share vers CGAL:

cd /usr/share/; ln -s /Users/Shared/CGAL-3.2.1 CGAL

Revenir ensuite aux sources de pgrouting. Puis faites

 ./configure --with-boost=/usr/local --with-cgal=/Users/Shared/CGAL-3.2.1 --with-gaul=/usr/local
make
sudo make install

En cas de problème avec le make

__Unwind_Resume
collect2: ld returned 1 exit status
make: *** [librouting.0.0.so] Error 

Éditer le makefile.in et remplacer la ligne

SHLIB_LINK += -lstdc++ $(TSP_LIBS) $(ALPHA_LIBS)

par

SHLIB_LINK += -lstdc++ $(TSP_LIBS) $(ALPHA_LIBS) -fexceptions

Ensuite

make clean
make
sudo make install 

A.4. Chargement des fonctionnalités de PgRouting

Pour créer une base routing, ayant les fonctionnalités de PgRouting, il nous suffira de faire

createdb routing
createlang plpgsql routing
psql -d routing -f [chemin_d_accces_vers]lwpostgis.sql
psql -d routing -f [chemin_d_accces_vers]spatial_ref_sys.sql
psql -d routing -f [chemin_d_accces_vers]routing.sql
psql -d routing -f [chemin_d_accces_vers]routing_postgis.sql

où [chemin_d_accces_vers] est le chemin d'accès vers le fichier attendu en fonction de votre OS (GNU/Linux ou Windows).

A.5. Importation d'un shapefile concernant des tronçons

Le shapefile que nous allons utiliser ici est fournit à l'adresse http://www.davidgis.fr/download/troncon_route.zip . Une fois le fichier dézippé, comme à l'accoutumée pour l'importer dans PostGIS,

shp2pgsql -DI troncon_route.shp troncon_route | psql routing

Les données géométriques du shapefile en question - visualisé dans QGIS- ressemblent à ca

Figure A.1. Le réseau avec 5 rond-points.

Le réseau avec 5 rond-points.

Pour le moment, ma table ne contient rien de bien intéressant. Dans la colonne sens, sens="sens direct" désignent un tronçon d'un rond-point.

SELECT gid,sens,astext(the_geom) FROM  troncon_route ORDER BY  gid "
 gid |    sens     |                astext
-----+-------------+--------------------------------------
   1 | double sens | MULTILINESTRING((1 0,5 0))
   2 | double sens | MULTILINESTRING((5 0,5 6))
   3 | double sens | MULTILINESTRING((0 7.5,3 7.5))
   4 | sens direct | MULTILINESTRING((3 7.5,3 7,4 6,5 6))
   5 | sens direct | MULTILINESTRING((5 6,6 6,7 7,7 7.5))
   6 | sens direct | MULTILINESTRING((7 7.5,7 8,6 9,5 9))
   7 | sens direct | MULTILINESTRING((5 9,4 9,3 8,3 7.5))
   8 | double sens | MULTILINESTRING((7 7.5,11 7.5))
   9 | double sens | MULTILINESTRING((11 7.5,11 11))
  10 | double sens | MULTILINESTRING((11 7.5,14 7.5))
  11 | double sens | MULTILINESTRING((14 7.5,21 7.5))
 ... | double sens | MULTILINESTRING((..................))

On aura compris que pour le moment, je ne dispose que des tronçons. Or pour la suite, il me faut définir les noeuds de mon réseau.

A.6. Obtention des noeuds du réseau

Les noeuds se présentent ainsi

Figure A.2. Les noeuds du réseau

Les noeuds du réseau

Ils ont été obtenus en utilisant les instructions SQL suivantes

BEGIN TRANSACTION;
-- ! ! !LA LIGNE SUIVANTE EST A DECOMMENTER SI ON DOIT RECHARGER CE BLOC D'INSTRUCTIONS SQL ! ! !
--SELECT drop_graph_tables('troncon_route');

/*
   Ajouter les colonnes adéquates
*/
ALTER TABLE troncon_route ADD column source_id int4;
ALTER TABLE troncon_route ADD column target_id int4;
ALTER TABLE troncon_route ADD column edge_id int4;
/*
   Mettre à jour le srid=1 sinon pgrouting gueule 8-(
*/
SELECT updategeometrysrid('troncon_route','the_geom',-1);

SELECT assign_vertex_id('troncon_route',0.00001);
/*
   Ok...Je crée mon graphe
*/
SELECT create_graph_tables('troncon_route', 'int4');

--SELECT UPDATE_cost_FROM_distance('troncon_route');
ALTER TABLE troncon_route_edges ADD column sens text;
ALTER TABLE troncon_route_edges ADD column x1 double precision;
ALTER TABLE troncon_route_edges ADD column y1 double precision;
ALTER TABLE troncon_route_edges ADD column x2 double precision;
ALTER TABLE troncon_route_edges ADD column y2 double precision;
ALTER TABLE troncon_route_edges ADD column edge_id int4;
/*
   Mise à jour des colonnes x1,y1,x2,y2 originaux par rapport aux données géométriques de la table troncon_route
   et mise à jour des colonnes sens et edge_id
*/
UPDATE troncon_route_edges SET cost=(SELECT length(the_geom) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom);
UPDATE troncon_route_edges SET x1=(SELECT x(startpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom);
UPDATE troncon_route_edges SET y1=(SELECT y(startpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom);
UPDATE troncon_route_edges SET x2=(SELECT x(endpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom);
UPDATE troncon_route_edges SET y2=(SELECT y(endpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom);
UPDATE troncon_route_edges SET edge_id=(SELECT edge_id FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.edge_id);
UPDATE troncon_route_edges SET sens=(SELECT sens::text FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.sens);
SELECT AddGeometryColumn( 'troncon_route_edges', 'the_geom', -1, 'MULTILINESTRING', 2 );
UPDATE troncon_route_edges SET the_geom=(SELECT the_geom FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom);
/*
    Tout ce qui est à double sens je le garde
*/
UPDATE troncon_route_edges SET reverse_cost=cost;
/*
    Paramétrer le coût des tronçons à sens unique
*/
UPDATE troncon_route_edges SET reverse_cost=-1 WHERE sens='sens direct';

END TRANSACTION;
VACUUM FULL ANALYZE ;

J'ai donc maintenant une table troncon_route_edges avec le contenu suivant.

SELECT id,sens,astext(the_geom),x1,y1,x2,y2,source,target,edge_id,cost,reverse_cost FROM  troncon_route_edges ORDER BY  id LIMIT 10

me renvoit

 id |    sens     |                astext                | x1 | y1  | x2 | y2  | source | target | edge_id |       cost       | reverse_cost
----+-------------+--------------------------------------+----+-----+----+-----+--------+--------+---------+------------------+--------------
  1 | double sens | MULTILINESTRING((1 0,5 0))           |  1 |   0 |  5 |   0 |      1 |      2 |       1 |                4 |            4
  2 | double sens | MULTILINESTRING((5 0,5 6))           |  5 |   0 |  5 |   6 |      2 |      3 |       2 |                6 |            6
  3 | double sens | MULTILINESTRING((0 7.5,3 7.5))       |  0 | 7.5 |  3 | 7.5 |      4 |      5 |       3 |                3 |            3
  4 | sens direct | MULTILINESTRING((3 7.5,3 7,4 6,5 6)) |  3 | 7.5 |  5 |   6 |      5 |      3 |       4 | 2.91421356237309 |           -1
  5 | sens direct | MULTILINESTRING((5 6,6 6,7 7,7 7.5)) |  5 |   6 |  7 | 7.5 |      3 |      6 |       5 | 2.91421356237309 |           -1
  6 | sens direct | MULTILINESTRING((7 7.5,7 8,6 9,5 9)) |  7 | 7.5 |  5 |   9 |      6 |      7 |       6 | 2.91421356237309 |           -1
  7 | sens direct | MULTILINESTRING((5 9,4 9,3 8,3 7.5)) |  5 |   9 |  3 | 7.5 |      7 |      5 |       7 | 2.91421356237309 |           -1
  8 | double sens | MULTILINESTRING((7 7.5,11 7.5))      |  7 | 7.5 | 11 | 7.5 |      6 |      8 |       8 |                4 |            4
  9 | double sens | MULTILINESTRING((11 7.5,11 11))      | 11 | 7.5 | 11 |  11 |      8 |      9 |       9 |              3.5 |          3.5
 10 | double sens | MULTILINESTRING((11 7.5,14 7.5))     | 11 | 7.5 | 14 | 7.5 |      8 |     10 |      10 |                3 |            3
(10 lignes)

A.7. Fonctionnalité shortest_path_astar()

Intéressons-nous en premier lieu à cette fonctionnalité du plus court chemin

A.7.1. Exemple pour les noeuds 38 et 48.

Je vais ici créer deux tables aller et retour

A.7.1.1. Pour l'aller (source=38 et target=48)

je fais

BEGIN TRANSACTION;

CREATE TABLE aller(gid int4) WITH oids;

SELECT AddGeometryColumn( 'aller', 'the_geom', -1, 'MULTILINESTRING', 2 );

INSERT INTO aller(the_geom) 
  (
   SELECT the_geom FROM troncon_route_edges  WHERE edge_id IN
     (SELECT edge_id FROM shortest_path_astar('SELECT  id,source,target,cost,
      reverse_cost, x1,y1,x2,y2 FROM troncon_route_edges',38,48,false,true)
     )
  );

END TRANSACTION;

ce qui me donnera le résultat suivant

Figure A.3. Parcours à l'aller

Parcours à l'aller

A.7.1.2. Pour le retour (source=48 et target=38)

Je fais

BEGIN TRANSACTION;

CREATE TABLE retour(gid int4) WITH oids;

SELECT AddGeometryColumn( 'retour', 'the_geom', -1, 'MULTILINESTRING', 2 );

INSERT INTO retour(the_geom) 
  (
   SELECT the_geom FROM troncon_route_edges  WHERE edge_id IN
     (SELECT edge_id FROM shortest_path_astar('SELECT  id,source,target,cost,
      reverse_cost, x1,y1,x2,y2 FROM troncon_route_edges',48,38,false,true)
     )
  );

END TRANSACTION;

Au niveau visualisation, on obtient

Figure A.4. Parcours au retour

Parcours au retour

A.7.2. Démo en ligne avec MapServer

Vous trouverez une démo en ligne de ce qui a été suggéré précédemment à l'adresse http://www.davidgis.fr/routing/

A.7.3. Tester soi-même la démo avec MapServer.

En pré-requis, il faut que vous ayez PhpMapScript et l'extension PostgreSQL pour Php d'activé. Vous trouverez un répertoire de ce tutoriel à http://www.davidgis.fr/download/pgrouting_demo.zip . Une fois le fichier téléchargé et dézippé, il suffit dans la mapfile mapfile/map.map d'adapter la portion de la partie WEB à votre configuration

  WEB
    IMAGEPATH "/var/www/tutorial/routing/tmp/" <--- A adapter
    IMAGEURL "/tutorial/routing/tmp/" <-- A adapter
      METADATA
      END
    QUERYFORMAT text/html
  END

et pour chaque layer, il faudra adapter la partie

"user=postgres dbname=routing password=empr888 host=localhost"

à votre propre configuration de PostgreSQL. Le shapefile troncon_route.shp est dans le sous-répertoire data

A.7.4. Tester sur un jeu de données réelles: jeu de tests GEOROUTE IGN

Sur le site de l'IGN, on peut avoir un échantillon des données gratuits. En remplissant un simple formulaire, on peut recevoir un CD contenant un jeu de test en zone agglomération sur la ville d'Orléans (45). Les données sont fournies au format MapInfo

A.7.4.1. Conversion du fichier MapInfo vers PostGIS

A la racine du CD, il faut se rendre à DONNEES/GEOROUTE/NAVIGATION/RESEAU_ROUTIER. Je vais utiliser ogr2ogr pour convertir une première fois le fichier en shapefile, puis en table PostGIS ensuite

cd /media/cdrom/DONNEES/GEOROUTE/NAVIGATION/RESEAU_ROUTIER
cp TRONCON_ROUTE.* ~/
cd
ogr2ogr -f "ESRI Shapefile" TRONCON_ROUTE.SHP TRONCON_ROUTE.TAB
shp2pgsql -dDI TRONCON_ROUTE.SHP troncon_route|iconv -f LATIN1 -t UTF-8|psql routing

A.7.4.2. Instructions SQL

Je vais effectuer les requêtes SQL suivantes

BEGIN TRANSACTION;

SELECT drop_graph_tables('troncon_route');

/*
   Ajouter les colonnes adéquates
*/
ALTER TABLE troncon_route ADD column source_id int4;
ALTER TABLE troncon_route ADD column target_id int4;
ALTER TABLE troncon_route ADD column edge_id int4;
/*
   Mettre à jour le srid=1 sinon pgrouting gueule 8-(
*/
SELECT updategeometrysrid('troncon_route','the_geom',-1);
UPDATE troncon_route SET the_geom=reverse(the_geom),sens='Sens direct' WHERE sens='Sens inverse';
SELECT assign_vertex_id('troncon_route',0.00001);
/*
   Ok...Je crée mon graphe
*/
SELECT create_graph_tables('troncon_route', 'int4');

--SELECT UPDATE_cost_FROM_distance('troncon_route');
ALTER TABLE troncon_route_edges ADD column sens text;
ALTER TABLE troncon_route_edges ADD column x1 double precision;
ALTER TABLE troncon_route_edges ADD column y1 double precision;
ALTER TABLE troncon_route_edges ADD column x2 double precision;
ALTER TABLE troncon_route_edges ADD column y2 double precision;
ALTER TABLE troncon_route_edges ADD column edge_id int4;
/*
   Mise à jour des colonnes x1,y1,x2,y2 originaux par rapport aux données géométriques de la table troncon_route
   et mise à jour des colonnes sens et edge_id
*/
SELECT update_cost_from_distance('troncon_route');

UPDATE troncon_route_edges SET x1=(SELECT x(startpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);

UPDATE troncon_route_edges SET y1=(SELECT y(startpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);

UPDATE troncon_route_edges SET x2=(SELECT x(endpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);

UPDATE troncon_route_edges SET y2=(SELECT y(endpoint(the_geom)) FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);

UPDATE troncon_route_edges SET edge_id=(SELECT edge_id FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.edge_id LIMIT 1);

UPDATE troncon_route_edges SET sens=(SELECT sens::text FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.sens LIMIT 1);

SELECT AddGeometryColumn( 'troncon_route_edges', 'the_geom', -1, 'MULTILINESTRING', 2 );

UPDATE troncon_route_edges SET the_geom=(SELECT the_geom FROM troncon_route g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);
/*
    Tout ce qui est à double sens je le garde
*/
UPDATE troncon_route_edges SET reverse_cost=cost;
/*
    Paramétrer le coût des tronçons à sens unique
*/
UPDATE troncon_route_edges SET reverse_cost=-1 WHERE sens='Sens direct';

END TRANSACTION;
VACUUM FULL ANALYZE ;

J'exporte ensuite la table troncon_route_edges en Shapefile en faisant

pgsql2shp -h localhost -u postgres -f troncon_route_edges.shp routing troncon_route_edges

A.7.4.3. Exemple

Depuis QGIS, je charge le shapefile troncon_route_edges.shp en deux fois: un pour avoir les sources (en fonction de x1,y1) et l'autre pour les target (en fonction de x2,y2). Pour celà, je joues sur les propriétés d'affichage que propose QGIS: clic-droit-> propriété, onglet "Etiquette" etc...Une fois tout celà fait, j'ai par exemple l'affichage suivant

Figure A.5. Une portion des noeuds de mon réseau (sources+targets)

Une portion des noeuds de mon réseau (sources+targets)

Je choisis ensuite un source et un target. Sur le CD, en me rendant à DONNEES/GEOROUTE_Raster/LZW/5K, je charge le fond gr_5k_ITagglo_v1.tif dans QGIS. Je crée ensuite les tables aller et retour comme fait précédemment. J'ai par exemple l'affichage suivant

Figure A.6. Parcours à l'aller

Parcours à l'aller

Figure A.7. Parcours au retour

Parcours au retour

A.8. Fonctionnalité TSP()

TSP signifie Traveling Salesman Problem, désignant ainsi le problème du voyageur du commerce. Pour une définition plus profane, il s'agit du problème du père Noë pour sa distribution de cadeauxl, ou de l'agent commercial qui doit visiter des clients dans un trajet professionnel en essayant d'optimiser son parcours - un peu de théorie des graphes, celà ne fait pas de mal que les puristes de la discipline m'excusent -. Le graphe que j'utilise ici est toujours le même. Nous conviendrons donc même si l'exemple qui sera fourni n'est pas des plus convaincants. Le but est en fade réussir à optimiser le parcours qui permet en passant par différents points attendus. Pour ma part ici, je vous renvois à la documentation réalisée par Gérald FENOY.

A.8.1. Exemple

Sur mon réseau, je prendrais source=7, comme point de départ. Je souhaite passer par les noeuds fournis ici dans un ordre quelconque 25,41,28,42.tsp() va me permettre de savoir quel est le meilleur ordre de parcours pour ces noeuds. J'éxecuerais donc la requête suivante

routing=# SELECT  * FROM tsp('SELECT  DISTINCT source AS source_id, x1 AS x, y1  AS y FROM
troncon_route_edges WHERE source IN (25,7,41,28,42)','25,7,41,28,42',7);

qui me renverra

 vertex_id | edge_id |         cost
-----------+---------+-----------------------
         7 |       8 | 5.76058619309876e-269
        41 |       8 | 5.76063240607483e-269
        28 |       8 |  5.7606786190509e-269
        42 |       8 | 6.83337886510596e-316
        25 |      32 | 2.13684976166112e-311
(5 lignes)

Les colonnes edge_id et cost n'ont pas d'importance. Seul compte la colonne vertex_id. Selon l'ordre énuméré des valeurs des lignes de la colonne de vertex_id, il me faudra donc pour optimiser mon parcours aller de 7 à 41, de 41 à 28, de 28 à 42 etc...

Et donc de proche en proche, il me suffira par exemple d'utiliser la fonctionnalité shortest_path_astar() pour reconstruire mon parcours ou d'utiliser la fonctionnalité tsp_astar_as_geometry_internal_id_directed(). Le parcours sera donc le suivant

Figure A.8. fonctionnalité TSP() pour les noeuds 7,41,28,42 et 25

fonctionnalité TSP() pour les noeuds 7,41,28,42 et 25

Pour le moment, nous allons nous limiter à la première solution proposée, à savoir utiliser shortest_path_astar() de proche en proche. Celà nou spermettra de voir les limites d'une telle solution. Il en découlera l'utilisation plus intéressante de tsp_astar_as_geometry_internal_id_directed()

A.8.2. Programme en C pour les appels successifs à shortest_path_astar()

En se basant sur l'exemple précédent la seule tâche qui nous resterait à effectuer serait de pouvoir automatiser les appels successifs à shortest_path_astar(). Pour ce faire, nous allons utiliser un programmetsp_trajet. Ce programme est écrit en C en se basant sur l'interface de programmation libpq de PostgreSQL. Il prendra comme paramètre la connexion à notre base, la liste des points à parcourir et le point de départ. L'utilisation de tsp_trajet sera par exemple

tsp_trajet "dbname=routing user=postgres"  "25,7,41,28,42,19" 7

Pour la simplifcation de l'exposé ici, le principe le plus simple sera de stocker dans une table trajet -que le programme se chargera d'alimenter - les divers parcours retournées par shortest_path_astar(). Notre table aura donc la structure suivante

#SELECT source,target,astext(the_geom) FROM trajet
 source | target |  the_geom
--------+--------+-----------------------
      7 |     41 | MULTILINESTRING((..))
     41 |     28 | MULTILINESTRING((..))
     28 |     42 | MULTILINESTRING((..))
     42 |     25 | MULTILINESTRING((..))
(4 lignes)

Par exemple, l'appel de shortest_path_astar() pour source=41 et target=28, correspondra à la requête

INSERT INTO trajet (source,target,the_geom) values 
(41,28,(SELECT geomunion(the_geom) FROM troncon_route_edges  WHERE edge_id IN 
     (SELECT edge_id FROM shortest_path_astar('SELECT  id,source,target,cost,reverse_cost,
     x1,y1,x2,y2 FROM troncon_route_edges',41,28,false,true))))

Le programme se chargera à chaque utilisation de recréer la table trajet. On notera ici l'emploi de GeomUnion() qui permet pour chaque couple

A.8.2.1. Le programme

On pourra proposer par exemple le programme suivant que l'on peut s'amuser à optimiser

#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <libpq-fe.h>

int main(int argc, char * argv[])
{
  PGresult *result,*result2;
  PGconn *conn;
     char query1[900]=
"CREATE OR REPLACE FUNCTION drop_geometrytable_if_exists(text) RETURNS text AS $$ \
DECLARE rec record;\
BEGIN \
IF NULLVALUE($1) THEN \
   RETURN 'ATTENTION: Table  non trouvée'; \
ELSE  \
     SELECT INTO rec tablename FROM pg_tables WHERE tablename = quote_ident($1);\
   IF FOUND THEN\
    EXECUTE 'SELECT DropGeometryTable('||quote_literal('public')||','||quote_literal($1)||')';\
    RETURN 'Effacement de la table ...OK';\
  END IF;\
 END IF;\
RETURN 'ATTENTION: Table non trouvée';END;$$\
 LANGUAGE plpgsql;\
SELECT drop_geometrytable_if_exists('trajet');\
CREATE TABLE trajet(gid int4,source int,target int) WITH oids;\
SELECT AddGeometryColumn( 'trajet', 'the_geom', -1, 'MULTILINESTRING', 2 );\
SELECT  vertex_id::int from tsp('select  distinct source as source_id\
, x1::double precision as x, y1::double precision as y from \
troncon_route_edges where source in (";
  
  if (argc != 4)
   {
    printf("Usage:%s \"dbname=XXX user=XXX host=XXX\ password=XXX\"  \"liste de points\" \
startpoint\n\nExemple %s \"dbname=routing user=postgres\"  \"43,30,18,41,1\" 30\n\n", argv[0],argv[0]);
       return EXIT_SUCCESS;
  }

  conn = PQconnectdb(argv[1]);

  if(PQstatus(conn) == CONNECTION_OK) {
    printf("Connexion OK\n");
    strcat(query1," ");
    strcat(query1,argv[2]);
    strcat(query1,")','");
   strcat(query1,argv[2]);
   strcat(query1,"',");
   strcat(query1,argv[3]);
   strcat(query1,")");
    //printf("%s\n",query1);

    result = PQexec(conn, query1);
    
   switch(PQresultStatus(result))
    {
      case PGRES_TUPLES_OK:
        {
             int r, n;
             int nrows = PQntuples(result);
             int nfields = PQnfields(result);
             printf("\n--------------------------------------------\n");
             printf(" %d insertions à effectuer dans la table projet",nrows); 
             printf("\n--------------------------------------------\n");
              for(r = 1; r < nrows; r++) 
              {
       for(n = 0; n < nfields; n++)
                {
           char	query2[800]="INSERT INTO trajet (source,target,the_geom) values (";
      strcat(query2,PQgetvalue(result, r-1, n));
                      strcat(query2,",");
      strcat(query2,PQgetvalue(result, r, n));
                      strcat(query2,",(");
 strcat(query2,"SELECT geomunion(the_geom) FROM troncon_route_edges\
  WHERE edge_id IN (SELECT edge_id FROM shortest_path_astar('SELECT  id,source,target,cost,reverse_cost,\
 x1,y1,x2,y2 FROM troncon_route_edges',");
      strcat(query2,PQgetvalue(result, r-1, n));
                      strcat(query2,",");
      strcat(query2,PQgetvalue(result, r, n));
                      strcat(query2,",false,true))))");
              //printf("%s\n",query2);
                      printf(" %s --> %s \n",PQgetvalue(result, r-1, n),PQgetvalue(result, r, n));
                      result2 = PQexec(conn, query2);
                }
               }
        }
    }
   PQclear(result);PQclear(result2);
  }
  else 
    printf("Echec de connexion: %s\n", PQerrorMessage(conn));

  PQfinish(conn);
  return EXIT_SUCCESS;
}

Pour la compilation, on pourra se baser sur le Makefile suivant à adapter selon vos besoins

INC=`pg_config --includedir`
LIB=`pg_config --libdir`

CFLAGS=-I$(INC)
LDLIBS=-L$(LIB) -lpq

ALL=tsp_trajet
all: $(ALL)

clean:
        @rm -f *.o *~ $(ALL)

A.8.3. Limites du programme

On l'aura compris de suite, ce genre de programme ne peut être utilisé dans un mileu de production et la gestion de la table trajet serait alors impossible dans un mode multi-utilisateur. Or PgRouting propose la fonctionnalité tsp_astar_as_geometry_internal_id_directed() qui permet de calculer à la volée les divers portions de parcours attendus.

A.9. Fonctionnalités shortest_path_astar_as_geometry_internal_id_directed() et tsp_astar_as_geometry_internal_id_directed()

Je profite ici de proposer deux exemples sur ces deux focntionnalités. Avant de commencer, nous allons importer un nouveau jeu de données issus de NavTeq. Nous en profiterons pour voir nos fonctionnalités graĉe à ce jeu.

A.9.1. Importation d'un jeu de données NavTeq

NavTeq on ne le présente plus. Tout le monde connaît. La société ADCI - http://www.adci.com - propose sur son site un jeu de données test sur le réseau routier à Washington DC. Pour se faire, il faut se rendre à http://www.adci.com/products/navteq/index.html. C'est surtout le jeu "NAVTEQ Premium - for routing applications" qu'il faut télécharger. Une simple inscription sur le site permet par la suite de pouvoir télécharger les données. Une fois téléchargées et décompressées, ce sera surtout le shapefile streets.shp que nous allons tester ici. Les données sont fournis dans le système GPS Mondial WGS 84

Le petit + c'est aussi la documentation qui accompagne les shapefiles. On y apprend comment NavTeq présente sa topologie, ses spécifités etc...Pour convertir le fichier en kml, rien de plus simple que

ogr2ogr -f KML ~/streets.kml streets.shp

Une petite importation du fichier streets.kml dans GoogleEarth nous renverra les screenshots suivantes

Figure A.9. GoogleEarth: réseau routier à Washington DC

GoogleEarth: réseau routier à Washington DC

ainsi que

Figure A.10. GoogleEarth: réseau routier à Washington DC

GoogleEarth: réseau routier à Washington DC

Pour l'importation des données, nous ferons tout simplement

shp2pgsql -DI streets.shp streets | psql routing

A.9.2. Noeuds du réseau et direction pour le routage

NavTeq propose déjà sa propre topologie. En vertu de la documentation qui accompagne les données, je pourrais bien me servir des champs ref_in_id et de nref_in_id pour les associer respectivement à mes champs source et target pour pgrouting. Or pour rester en conformité avec ma documentation, je préfère reconstruit les noeuds. De plus pour le sens de parcours des tronçons, il faut savoir que ma table streets contient le champ dir_travel (Direction Travel) qui a les spécifités suivantes

  1. dir_travel='B' c'est une tronçon à double sens;

  2. dir_travel='F' c'est une tronçon à sens unique. Le sens de parcoursde l'arc est le sens direct par rapport au noeud de référence (le point de départ est le premier point de mon arc);

  3. dir_travel='T' c'est une tronçon à sens unique. Le sens de parcours de l'arc est le sens opposé par rapport au noeud de référence (le point de départ est le dernier point de mon arc);

Pour la suite, je vais donc appliquer le même principe que celui des instructions SQL que j'ai employé pour le jeu de données GEOROUTE de l'IGN. Par analogie, la colonne dir_travel jouera donc le même rôle que celle de la colonne sens du GEOROUTE: Les différences notables seront donc

BEGIN TRANSACTION;
ALTER TABLE streets ADD column source_id int4;
ALTER TABLE streets ADD column target_id int4;
ALTER TABLE streets ADD column edge_id int4;
SELECT updategeometrysrid('streets','the_geom',-1);
UPDATE streets SET the_geom=reverse(the_geom),dir_travel='F' WHERE dir_travel='T';
SELECT assign_vertex_id('streets',0.00001);
SELECT create_graph_tables('streets', 'int4');
ALTER TABLE streets_edges ADD column dir_travel text;
ALTER TABLE streets_edges ADD column x1 double precision;
ALTER TABLE streets_edges ADD column y1 double precision;
ALTER TABLE streets_edges ADD column x2 double precision;
ALTER TABLE streets_edges ADD column y2 double precision;
ALTER TABLE streets_edges ADD column edge_id int4;
SELECT update_cost_from_distance('streets');
UPDATE streets_edges SET x1=(SELECT x(startpoint(the_geom)) FROM streets g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);
UPDATE streets_edges SET y1=(SELECT y(startpoint(the_geom)) FROM streets g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);
UPDATE streets_edges SET x2=(SELECT x(endpoint(the_geom)) FROM streets g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);
UPDATE streets_edges SET y2=(SELECT y(endpoint(the_geom)) FROM streets g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);
UPDATE streets_edges SET edge_id=(SELECT edge_id FROM streets g WHERE g.edge_id=id GROUP BY id,g.edge_id LIMIT 1);
UPDATE streets_edges SET dir_travel=(SELECT dir_travel::text FROM streets g WHERE g.edge_id=id GROUP BY id,g.dir_travel LIMIT 1);
SELECT AddGeometryColumn( 'streets_edges', 'the_geom', -1, 'MULTILINESTRING', 2 );
UPDATE streets_edges SET the_geom=(SELECT the_geom FROM streets g WHERE g.edge_id=id GROUP BY id,g.the_geom LIMIT 1);
UPDATE streets_edges SET reverse_cost=cost;
UPDATE streets_edges SET reverse_cost=-1 WHERE dir_travel='F';
END TRANSACTION;

A.9.3. Modifications nécessaires sur la table streets_edges

Il faut effectuer les commandes suivantes

ALTER TABLE street_edges RENAME id TO gid;
ALTER TABLE street_edges DROP COLUMN edge_id;
alter TABLE street_edges RENAME cost TO length;

Ces modifications sont nécessaires notamment pour l'emploi des fonctionnalités. On en profitera aussi pour créer les index sur les colonnes source et target

CREATE INDEX k1 on streets_edges(source);
CREATE INDEX k2 on streets_edges(target);
VACUUM FULL ANALYZE;

A.9.4. Exemple avec shortest_path_astar_as_geometry_internal_id_directed()

Comme mon graphe est orienté, pour aller par exemple du noeud source=5274 au noeud target=5488, il me faudra faire

SELECT * FROM shortest_path_astar_as_geometry_internal_id_directed('street_edges',5274,5488,false,true);

qui renverra les champs gid et the_geom. Ce qui implique qu'elle puisse être directement utilisée par MapServer. Sur les deux images suivantes, les tronçons à double sens (respectivement à sens unique) sont en bleu (respectivement en vert).

Figure A.11. Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5274 et target=5488 (sens aller)

Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5274 et target=5488 (sens aller)

Figure A.12. Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5274 et target=5488 (sens aller)

Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5274 et target=5488 (sens aller)

et pour le retour

SELECT * FROM shortest_path_astar_as_geometry_internal_id_directed('street_edges',5488,5274,false,true);

qui renverra comme résultat

Figure A.13. Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5488 et target=5274 (sens retour)

Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5488 et target=5274 (sens retour)

Figure A.14. Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5488 et target=5274 (sens retour)

Fonction shortest_path_astar_as_geometry_internal_id_directed() sur source=5488 et target=5274 (sens retour)

A.9.5. Exemple avec tsp_astar_as_geometry_internal_id_directed()

Comme mon graphe est orienté, je vais ici considérer les noeuds 5403,5822,338,7106,6043,1952. Je prendrais comme point de départ source=1952

Nous obtiendrons directement les enregistrements the_geom en faisant

SELECT the_geom FROM tsp_astar_as_geometry_internal_id_directed('streets_edges','5403,5822,338,7106,6043,1952',1952,0.03,false,true);

qui renverra les champs gid et the_geom. Ce qui implique qu'elle puisse être directement utilisée par MapServer.

Figure A.15. Fonction tsp_astar_as_geometry_internal_id_directed()

Fonction tsp_astar_as_geometry_internal_id_directed()

Une exportation de cette requête dans un fichier kml par

ogr2ogr -f KML ~/parcours.kml PG:'host=localhost dbname=routing user=postgres' \
-sql "SELECT (dump).geom FROM (SELECT dump(the_geom) FROM \
tsp_astar_as_geometry_internal_id_directed('streets_edges','5403,5822,338,7106,6043,1952',1952,.03,false,true)) AS foo"

me renverra l'image suivante dans GoogleEarth

Figure A.16. Fonction tsp_astar_as_geometry_internal_id_directed() vue depuis GoogleEarth

Fonction tsp_astar_as_geometry_internal_id_directed() vue depuis GoogleEarth