L'algorythme de Dijkstra est connu en informatique (théorie des graphes) pour permettre de trouver le parcours du meilleur chemin. Porté initialement comme une extension dans le projet CARTOWEB http://www.cartoweb.org sous le nom de pgdijkstra, il permet en autre de trouver sur un réseau routier par exemple de trouver le meilleur chemin à parcourir pour se rendre d'un point à un autre

Ici je vais me permettre de partir des données fournies dans les tables small_roads, et great_roads du chapitre "Tutoriaux", de tronçonner les rues et de les mettre dans une tables roads dont voici le contenu

BEGIN;
CREATE TABLE "roads" (gid serial, "source_id" int4, "target_id" int4,"orientation" varchar);
SELECT AddGeometryColumn('','roads','the_geom',-1,'MULTILINESTRING',2);
--
-- Rue Figaro en tronçons
--
 INSERT INTO roads values (0,1,2,'double sens',GeometryFromText('MULTILINESTRING((60 87.5,60 115))',-1));
 INSERT INTO roads values (1,2,3,'double sens',GeometryFromText('MULTILINESTRING((60 115,85 120))',-1));
--
-- Rue Voltaire
--
 INSERT INTO roads values (2,1,4,'sens inverse',GeometryFromText('MULTILINESTRING((60 87.5,60 25))',-1));
--
-- Rue Paul Valéry en tronçons
-- 
 INSERT INTO roads values (3,5,6,'double sens',GeometryFromText('MULTILINESTRING((1 1,20 23))',-1));
 INSERT INTO roads values (4,6,4,'double sens',GeometryFromText('MULTILINESTRING((20 23,60 25))',-1));
 INSERT INTO roads values (5,4,7,'double sens',GeometryFromText('MULTILINESTRING((60 25,85 36))',-1));
--
-- Rue du Général de Gaulle en tronçons
-- 
 INSERT INTO roads values (6,8,1,'double sens',GeometryFromText('MULTILINESTRING((1 87.5,60 87.5))',-1));
 INSERT INTO roads values (7,1,9,'double sens',GeometryFromText('MULTILINESTRING((60 87.5,85 87.5))',-1));
 INSERT INTO roads values (8,9,10,'double sens',GeometryFromText('MULTILINESTRING((85 87.5,150 87.5))',-1));
--
-- Rue Aristide Briand en tronçons
-- 
 INSERT INTO roads values (9,11,3,'double sens',GeometryFromText('MULTILINESTRING((85 135,85 120))',-1));
 INSERT INTO roads values (10,3,9,'double sens',GeometryFromText('MULTILINESTRING((85 120,85 87.5))',-1));
 INSERT INTO roads values (11,9,7,'double sens',GeometryFromText('MULTILINESTRING((85 87.5,85 36))',-1));
 INSERT INTO roads values (12,7,12,'double sens',GeometryFromText('MULTILINESTRING((85 36,85 1))',-1));

ALTER TABLE ONLY "roads" ADD CONSTRAINT "roads_pkey" PRIMARY KEY (gid);
ALTER TABLE roads ADD COLUMN edge_id integer;
END;

Voici le contenu de ma table

SELECT gid,source_id,target_id,orientation,astext(the_geom) from roads;

 gid | source_id | target_id | orientation  |               astext
-----+-----------+-----------+--------------+-------------------------------------
   0 |         1 |         2 | double sens  | MULTILINESTRING((60 87.5,60 115))
   1 |         2 |         3 | double sens  | MULTILINESTRING((60 115,85 120))
   2 |         1 |         4 | sens inverse | MULTILINESTRING((60 87.5,60 25))
   3 |         5 |         6 | double sens  | MULTILINESTRING((1 1,20 23))
   4 |         6 |         4 | double sens  | MULTILINESTRING((20 23,60 25))
   5 |         4 |         7 | double sens  | MULTILINESTRING((60 25,85 36))
   6 |         8 |         1 | double sens  | MULTILINESTRING((1 87.5,60 87.5))
   7 |         1 |         9 | double sens  | MULTILINESTRING((60 87.5,85 87.5))
   8 |         9 |        10 | double sens  | MULTILINESTRING((85 87.5,150 87.5))
   9 |        11 |         3 | double sens  | MULTILINESTRING((85 135,85 120))
  10 |         3 |         9 | double sens  | MULTILINESTRING((85 120,85 87.5))
  11 |         9 |         7 | double sens  | MULTILINESTRING((85 87.5,85 36))
  12 |         7 |        12 | double sens  | MULTILINESTRING((85 36,85 1))
(13 lignes)

On va maintenant utiliser la fonction assign_vertex_id() en lui appliquant la distance de 0.01

SELECT Assign_Vertex_Id('roads',0.01);

On crée maintenant notre graphe grâce à la fonction create_graph_tables():

SELECT Create_Graph_Tables('roads','int4');

qui créera aussitôt deux tables supplémentaires roads_vertices et roads_edges dont les contenus respectifs sont

testgis=# SELECT * from roads_vertices;
 id | geom_id
----+---------
  1 |       1
  2 |       2
  3 |       3
  4 |       4
  5 |       5
  6 |       6
  7 |       7
  8 |       8
  9 |       9
 10 |      10
 11 |      11
 12 |      12
(12 lignes)

et

testgis=# SELECT * from roads_edges;
 id | source | target | cost | reverse_cost
----+--------+--------+------+--------------
  1 |      1 |      2 |      |
  2 |      2 |      3 |      |
  3 |      1 |      4 |      |
  4 |      5 |      6 |      |
  5 |      6 |      4 |      |
  6 |      4 |      7 |      |
  7 |      8 |      1 |      |
  8 |      1 |      9 |      |
  9 |      9 |     10 |      |
 10 |     11 |      3 |      |
 11 |      3 |      9 |      |
 12 |      9 |      7 |      |
 13 |      7 |     12 |      |
(13 lignes)

Pour l'instant, cette dernière table a sa colonne des coûts (colone cost) et aussi des coûts inversés (colonne reverse_cost) vides. Nous allons les remplir.

Figure A.2. Graphe orienté, la rue Voltaire est à sens unique (en jaune sur le graphe,du noeud 4 au noeud 1)

Nous allons baser les coûts pour se rendre d'un noeud à un autre (voir figure plus loin dans le document) sur les distances qui les séparent. L'équipe de CartoWeb a ainsi prévu à cette fin, la fonction update_cost_from_distance() - repris par l'équipe de PgRouting -

SELECT update_cost_from_distance('roads');

Nous allons aussi estimé que notre graphe est orienté et que l'on peut aussi pour certains points hormis le 1 et le 4 aller d'un point à un autre.

UPDATE roads_edges SET reverse_cost=cost;

ce qui renverra donc

testgis=# SELECT *  from roads_edges order by id;
 id | source | target |       cost       |   reverse_cost
----+--------+--------+------------------+------------------
  1 |      1 |      2 |             27.5 |             27.5
  2 |      2 |      3 | 25.4950975679639 | 25.4950975679639
  3 |      1 |      4 |             62.5 |             62.5
  4 |      5 |      6 | 29.0688837074973 | 29.0688837074973
  5 |      6 |      4 | 40.0499687890016 | 40.0499687890016
  6 |      4 |      7 | 27.3130005674953 | 27.3130005674953
  7 |      8 |      1 |               59 |               59
  8 |      1 |      9 |               25 |               25
  9 |      9 |     10 |               65 |               65
 10 |     11 |      3 |               15 |               15
 11 |      3 |      9 |             32.5 |             32.5
 12 |      9 |      7 |             51.5 |             51.5
 13 |      7 |     12 |               35 |               35
(13 lignes)

Pour que la Rue Voltaire soit à sens unique ( de 4 --> 1), il faut lui appliquer un coût relativement élévé => cost = 1000000: 4 ->1 sens inverse.

D'où la commande

UPDATE roads_edges SET cost=1000000 WHERE id=3

Ce qui donne pour la table roads_edges:

# SELECT * from roads_edges ORDER BY 1;
 id | source | target |       cost       |   reverse_cost
----+--------+--------+------------------+------------------
  1 |      1 |      2 |             27.5 |             27.5
  2 |      2 |      3 | 25.4950975679639 | 25.4950975679639
  3 |      1 |      4 |          1000000 |             62.5
  4 |      5 |      6 | 29.0688837074973 | 29.0688837074973
  5 |      6 |      4 | 40.0499687890016 | 40.0499687890016
  6 |      4 |      7 | 27.3130005674953 | 27.3130005674953
  7 |      8 |      1 |               59 |               59
  8 |      1 |      9 |               25 |               25
  9 |      9 |     10 |               65 |               65
 10 |     11 |      3 |               15 |               15
 11 |      3 |      9 |             32.5 |             32.5
 12 |      9 |      7 |             51.5 |             51.5
 13 |      7 |     12 |               35 |               35
(13 lignes)

Nous utiliserons ici la fonction shortest_path(). Nous partirons du principe que le graphe est orienté et que nous utilisons aussi les couts inversés.

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.7. Parcours à l'aller

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.8. Parcours au retour

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

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=1000000 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

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.9. 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.10. Parcours à l'aller

Figure A.11. Parcours au retour

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)