C’est un sujet que j’avais déjà traité dans un post antérieur, mais comme le blog de l’époque a disparu, il a fallu que je me replonge dans les bas-fonds de postgresql, postgis et pgRouting. Heureusement que les articles de Anita Graser (alias underdarck) sont encore en ligne, eux, mais ils ont nécessité quelques ajustements avec pgRouting 2!

Si vous en êtes là c’est que vous cherchez un moyen de faire des calculs d’isochrones et que, petit à petit, la solution pgRouting (même si elle peut sembler disproportionnée) s’impose à vous.

A ce stade, je partirais du principe que vous avez un serveur postgreSQL/postGIS fonctionnel. Vous trouverez plusieurs ressources sur Internet. Pour ma part, j’ai réalisé ce petit tuto sur Fedora 20 (un bon tuto pour l’installation postgresql postgis su le blog de Marcus Neteler et pgRouting est dans les packages de la distribution).

Attention toutefois à la version. Il se trouve que, sur Fedora, pgRouting est encore à sa version 1.05 qui nécessite quelques ajustement à la mano (révélé par Gregor the map guy). De mon côté, j’ai plutôt opté pour une compilation depuis les sources qui sont disponibles sur gitHub.

Créer un réseau « routable »

La première chose à faire est de préparer la base de données :

createdb --encoding=UTF8 --owner=delaye geodb
psql -d geodb -c "CREATE EXTENSION postgis;"
psql -d geodb -c "CREATE EXTENSION postgis_topology;"
psql -d geodb -c "CREATE EXTENSION pgrouting;"

Ensuite,  il faut importer un shapefile dans postgreSQL. Vous pouvez utiliser un client graphique comme shp2pgsql-gui ou deux petites lignes de commande (des infos ici )

#transformer en requête SQL un shapefile
shp2pgsql -I -t 2D -s 2154 reseau_total_AOC.shp acces_banyuls access_banyuls > acces_banyuls.sql
#envoyer les données dans la base 
psql -f acces_banyuls.sql -d geodb

Vous pouvez dès à present vous connecter à la base avec QGIS et avoir accès aux données

Pour pouvoir créer un graph connecté, il faut identifier pour chaque ligne ou polyligne le début et la fin. C’est ce que nous allons faire avec :

--evaluer  le type de la colonne geom
select * from geometry_columns;
--pour évaluer  la dimension des coordonnées de vos géométries
select distinct st_coordDim(geom) from acces_banyuls;

--convertir la colonne en deux dimensions si les données sont en 2D si vous n'avez pas utilisé l'option -t 2D de shp2pgsql
alter table acces_banyuls alter column geom type geometry(Multilinestring, 2154) using st_force2D(geom);

On va ensuite devoir transformer des multiLine-string en line-string dans postGIS parce que dans sa version 2 pgRouting ne travaille plus qu’avec ça .

--d'utilisation de st_dump pour transformer multiLine-string en line-string
CREATE TABLE simple_ways AS SELECT *, (ST_Dump(geom)).geom AS simple_geom FROM acces_banyuls WHERE geom IS NOT NULL;
--au passage on suprime la colonne multiline string
ALTER TABLE simple_ways DROP COLUMN geom;

On va ensuite créer une colonne temps de trajet à partir d’une colonne existante dans la base de données qui contient, pour chaque tronçon, la vitesse limite :

--mise à jour des temps de trajet
ALTER TABLE simple_ways ADD COLUMN length_m INTEGER;
UPDATE simple_ways SET length_m = st_length(st_transform(simple_geom,2154));

ALTER TABLE simple_ways ADD COLUMN traveltime_min DOUBLE PRECISION;
UPDATE simple_ways set traveltime_min = length_m / (vtss_pt / 60 * 1000)  --pour avoir la vitesse en metre

On peut ensuite essayer les fonctions de pgRouting proposer par Anita Graser

-- ajouter des col "source" et "target" à la table
ALTER TABLE acces_banyuls ADD COLUMN "source" integer;
ALTER TABLE acces_banyuls ADD COLUMN "target" integer;

--Creation de colonne startPoint et endpoint (syntaxe à changer ne pas oublier st_)
CREATE OR REPLACE VIEW road_ext AS
   SELECT *, st_startpoint(simple_geom), st_endpoint(simple_geom)
   FROM simple_ways;

--Une table qui contient tout les noeuds
CREATE TABLE node AS
   SELECT row_number() OVER (ORDER BY foo.p)::integer AS id,
          foo.p AS the_geom
   FROM (        
      SELECT DISTINCT road_ext.st_startpoint AS p FROM road_ext
      UNION
      SELECT DISTINCT road_ext.st_endpoint AS p FROM road_ext
   ) foo
   GROUP BY foo.p;

   
--creation du réseau  -> 128330 ms 
CREATE TABLE network AS
   SELECT a.*, b.id as start_id, c.id as end_id
   FROM road_ext AS a
      JOIN node AS b ON a.st_startpoint = b.the_geom
      JOIN node AS c ON a.st_endpoint = c.the_geom;

Pour la dernière commande il faut faire attention, en effet à chaque fois que vous regénérez une table node, les id de ces derniers peuvent changer. Donc quand vous allez devoir spécifier l’ID du noeud source (6791 dans la requête suivante), faites attention de bien le mettre à jour.

--Creation de la table qui contient les temps de trajet
--besoin d'une colonne avec travelTime pour faire le calcul
-- temps de calcul 123190 ms
create table cost_road as
select
    id,
    the_geom,
    (select sum(cost) from (
       SELECT * FROM pgr_dijkstra('
       SELECT gid AS id,
          start_id::int4 AS source,
          end_id::int4 AS target,
          traveltime_min::float8 AS cost
       FROM network',
       6791,
       id,
       false,
       false)) as foo ) as cost
from node;

Vous pouvez ensuite charger ces points dans Qgis et lancer une interpolation avec l’extension interpolation, par exemple, qui va vous en proposer deux types :

1. TIN (Triangulated Irregular Netwaork)
2. IDW (Inverse Distance Weighted)

Voilà le résultat d’une interpolation IDW (avec une discrétisation moche je vous l’accorde 🙂 mais on a qu’à dire que c’est pas l’objet ici).

Ou, si vous voulez une interpolation plus propre, vous pouvez aussi passer sous R.