[PostGIS/PgRouting] With french routable data

Road-network isochrones with PostGIS & pgRouting (Docker)

This post documents a complete workflow to compute road-network isochrones using PostGIS and pgRouting, running inside Docker.
It used Data (BDTopo Dataset) from the IGN, France’s national mapping agency.
It simulates the time taken by fire trucks to reach different locations from a set of origin points (fire stations).

⚠️ Important disclaimer (read first or don't)

• speed information in BD TOPO is not always correct, and optimistic, work is being done to make it more accurate.
• intersections are not always modelled correctly (see pgr_separateTouching and/or use pgr_separateCrossing)
  • https://docs.pgrouting.org/latest/en/pgRouting-concepts.html#id65
• snapping is approximate (For the starting point, we get the closest point on the graph, it can be far away, or not truly be accessible in reality)
• roads are sometimes included entirely even if only partially reachable

The result should be understood as a network reach envelope, not a real travel-time model.

1. PostgreSQL + pgRouting in Docker

docker run -d \
  --name pgrouting \
  -p 5432:5432 \
  -e POSTGRES_DB=gis \
  -e POSTGRES_USER=gis \
  -e POSTGRES_PASSWORD=gis \
  -v $(pwd)/pgdata:/var/lib/postgresql \
  pgrouting/pgrouting:18-3.6-4.0

CREATE EXTENSION postgis;
CREATE EXTENSION pgrouting;

2. Preparing the Road Network

Upload road network to the database (gdal, qgis). The data used here was downloaded from the french IGN https://geoservices.ign.fr/bdtopo

ST_Force2D is used to drop Z/M values if present.

CREATE TABLE bdtopo.edges AS
SELECT
  id::bigint AS id,
  ST_Force2D(geom)::geometry(LineString, 2154) AS geom,
  vitesse_moyenne_vl
FROM bdtopo.troncon_de_route
WHERE geom IS NOT NULL;

Always add an index on the geometry column for performance.

CREATE INDEX edges_gix ON bdtopo.edges USING GIST (geom);

Add required columns for pgRouting.

ALTER TABLE bdtopo.edges
  ADD COLUMN source bigint,
  ADD COLUMN target bigint,
  ADD COLUMN cost double precision,
  ADD COLUMN reverse_cost double precision;

3. Travel Time (Length × Speed)

Replace NULL speeds with 0 (unusable edges).

UPDATE bdtopo.edges
SET vitesse_moyenne_vl = 0
WHERE vitesse_moyenne_vl IS NULL;

Compute costs based on length and speed (in km/h).
With BD TOPO, I exaggerated the cost by 15% to get more optimistic results.
Nullif is used to avoid division by zero; 1000.0 / 3600.0 converts km/h to m/s, I hope

UPDATE bdtopo.edges
SET cost = 1.15 *ST_Length(geom)
           / NULLIF(vitesse_moyenne_vl::double precision * 1000.0 / 3600.0, 0);

UPDATE bdtopo.edges
SET reverse_cost = 1.15 * ST_Length(geom)
           / NULLIF(vitesse_moyenne_vl::double precision * 1000.0 / 3600.0, 0);

Costs are expressed in seconds.

Intersections are free (if a node exists ofc), speeds are approximate, results are optimistic.

4. Graph Vertices

This is the point layer representing all start/end points of edges.

CREATE TABLE bdtopo.edges_vertices AS
SELECT *
FROM pgr_extractVertices(
  $$SELECT id, geom FROM bdtopo.edges ORDER BY id$$
);

Always add indexes for performance, well not always but when it makes sense.

CREATE INDEX edges_vertices_gix ON bdtopo.edges_vertices USING GIST (geom);
CREATE INDEX edges_vertices_id  ON bdtopo.edges_vertices (id);

For each edge (line), find the start and end point and link to the corresponding vertex id (calculated previously).

UPDATE bdtopo.edges e
SET source = v.id
FROM bdtopo.edges_vertices v
WHERE ST_StartPoint(e.geom) = v.geom;

UPDATE bdtopo.edges e
SET target = v.id
FROM bdtopo.edges_vertices v
WHERE ST_EndPoint(e.geom) = v.geom;

The result is visible in QGIS.

cost and reverse_cost are in seconds.

5. Snapping Origin Points

In the example we use fire stations as origin points (these are not the real locations).
You only need a table with points and an id column.

ALTER TABLE bdtopo.fire_stations
  ADD COLUMN IF NOT EXISTS start_vid bigint,
  ADD COLUMN IF NOT EXISTS snap_dist_m double precision;

We need to snap them to the closest vertex in the road network.
This is an approximation; ideally we would snap to the closest point on the closest edge.

UPDATE bdtopo.fire_stations p
SET
  start_vid = (
    SELECT v.id
    FROM bdtopo.edges_vertices v
    ORDER BY v.geom <-> p.geom
    LIMIT 1
  ),
  snap_dist_m = (
    SELECT ST_Distance(p.geom, v.geom)
    FROM bdtopo.edges_vertices v
    ORDER BY v.geom <-> p.geom
    LIMIT 1
  );

As we can see, the closest vertex is not very representative of the actual road network location.

6. Isochrone Computation

The aim here is to compute polygons representing areas reachable within a set of time thresholds (5, 10, 15, 20 minutes).

First create the output table.

CREATE TABLE bdtopo.isochrones (
  fire_stations_id bigint,
  minutes int,
  geom geometry(Polygon, 2154)
);

And here is the complete query to compute the isochrones.

TRUNCATE bdtopo.isochrones;

WITH dd20 AS (
    SELECT
        p.id AS fire_station_id,
        d.node,
        d.agg_cost
    FROM bdtopo.sis p
     JOIN LATERAL pgr_drivingDistance(
        $$SELECT id, source, target, cost, reverse_cost FROM bdtopo.edges$$,
        p.start_vid,
        20 * 60,
                                              directed := false
        ) d ON TRUE
WHERE p.start_vid IS NOT NULL -- and p.id = 7 TEST ONLY
  AND d.agg_cost <= 20*60
    ),
    bucketed AS (
SELECT
    fire_station_id,
    node,
    CASE
    WHEN agg_cost <  5*60 THEN 5
    WHEN agg_cost < 10*60 THEN 10
    WHEN agg_cost < 15*60 THEN 15
    ELSE 20
    END AS minutes
FROM dd20
    ),
    reach_pts AS (
SELECT
    b.fire_station_id,
    b.minutes,
    v.geom
FROM bucketed b
    JOIN bdtopo.edges_vertices v ON v.id = b.node
    )
INSERT INTO bdtopo.isochrones (fire_station_id, minutes, geom)
SELECT
    fire_station_id,
    minutes,
    ST_Buffer(
            ST_ConcaveHull(ST_Collect(geom), 0.10),
            5
    )::geometry(Polygon,2154) AS geom
FROM reach_pts
GROUP BY fire_station_id, minutes;

Quick tip: As polygons are not donuts.

You can still use QGIS to make sure the minutes taken are visible:

• click on layer properties → Symbology → Categorized
• check Control feature rendering order

Here is what it looks like in QGIS :)

And vs OSM isochrone tool (ORS tools), which handles a lot more (turns, roundabouts, etc.)

For 5 and 10 minutes both sets of data are quite close, but this gets worse the further the distance.

Some explanation about the above query

For each point in bdtopo.fire_stations, compute travel-time isochrones for:

• 5 minutes
• 10 minutes
• 15 minutes
• 20 minutes

Using a road network stored in bdtopo.edges.

Costs are expressed in seconds.

SQL Workflow

-- TRUNCATE bdtopo.isochrones;

WITH dd20 AS (
    SELECT
        p.id AS fire_station_id,
        d.node,
        d.agg_cost
    FROM bdtopo.sis p
     JOIN LATERAL pgr_drivingDistance(
        $$SELECT id, source, target, cost, reverse_cost FROM bdtopo.edges$$,
        p.start_vid,
        20 * 60,
        directed := false
        ) d ON TRUE
    WHERE p.start_vid IS NOT NULL
      AND d.agg_cost <= 20*60
),
bucketed AS (
    SELECT
        fire_station_id,
        node,
        CASE
            WHEN agg_cost <  5*60 THEN 5
            WHEN agg_cost < 10*60 THEN 10
            WHEN agg_cost < 15*60 THEN 15
            ELSE 20
        END AS minutes
    FROM dd20
),
reach_pts AS (
    SELECT
        b.fire_station_id,
        b.minutes,
        v.geom
    FROM bucketed b
    JOIN bdtopo.edges_vertices v ON v.id = b.node
)
INSERT INTO bdtopo.isochrones (fire_station_id, minutes, geom)
SELECT
    fire_station_id,
    minutes,
    ST_Buffer(
        ST_ConcaveHull(ST_Collect(geom), 0.10),
        5
    )::geometry(Polygon,2154) AS geom
FROM reach_pts
GROUP BY fire_station_id, minutes;

Step-by-Step Explanation

1. Compute reachable vertices within 20 minutes (dd20)

• pgr_drivingDistance is executed once per SIS point
• Maximum distance: 20 minutes (1200 seconds)
• Output nodes are filtered to agg_cost <= 20*60

Result:

(fire_station_id, node_id, travel_time_seconds)

2. Assign vertices to time bands (bucketed)

Vertices are classified into exclusive time ranges:

Travel time	Bucket
0–4:59	5 min
5–9:59	10 min
10–14:59	15 min
15–20:00	20 min

3. Convert vertices to points (reach_pts)

Each reachable node is joined to bdtopo.edges_vertices
to retrieve its POINT geometry.

4. Build polygons from points

For each (fire_station_id, minutes) group:

1. Points are collected
2. A tight concave hull is computed using alpha 0.10
3. A small buffer (5 m) smooths the geometry

Interpretation

These polygons represent:

Approximate envelopes of reachable network vertices

They are suitable for:

• exploratory analysis
• accessibility screening
• visual communication

They are not precise travel-time surfaces.

Tip

If you want cumulative isochrones (0–5, 0–10, 0–15, 0–20),
And make sur the smaller isochrones are drawn on top of the larger ones.

7. Visualization with points

Here is a small snippet to extract isochrone points for visualization,
and see for each point which category (5, 10, 15, 20 minutes) it belongs to.

DROP TABLE IF EXISTS bdtopo.isochrone_points;
CREATE TABLE bdtopo.isochrone_points (
  fire_station_id   bigint,
  minutes  int,
  vid      bigint,
  agg_cost double precision,
  geom     geometry(Point, 2154)
);

CREATE INDEX isochrone_points_gix ON bdtopo.isochrone_points USING GIST (geom);
CREATE INDEX isochrone_points_idx ON bdtopo.isochrone_points (fire_station_id, minutes);

INSERT INTO bdtopo.isochrone_points (fire_station_id, minutes, vid, agg_cost, geom)
WITH dd20 AS (
  SELECT
    p.id AS fire_station_id,
    d.node,
    d.agg_cost
  FROM bdtopo.sis p
  JOIN LATERAL pgr_drivingDistance(
    $$SELECT id, source, target, cost, reverse_cost FROM bdtopo.edges$$,
    p.start_vid,
    20 * 60,
    directed := false
  ) d ON TRUE
  WHERE p.start_vid IS NOT NULL -- and p.id = 7 TEST ONLY
),
bucketed AS (
  SELECT
    fire_station_id,
    node,
    agg_cost,
    CASE
      WHEN agg_cost <=  5 * 60 THEN  5
      WHEN agg_cost <= 10 * 60 THEN 10
      WHEN agg_cost <= 15 * 60 THEN 15
      ELSE 20
    END AS minutes
  FROM dd20
)
SELECT
  b.fire_station_id,
  b.minutes,
  b.node AS vid,
  b.agg_cost,
  v.geom
FROM bucketed b
JOIN bdtopo.edges_vertices v
  ON v.id = b.node;

You can visualize the result in QGIS.
This helps you understand what can be wrong/explained with/by the network data.

Conclusion

These isochrones are approximations.
They are useful for exploration, not for precision routing—unless you have a very high-quality road network. :)

I did not work on the base data itself, it would need to be preprocessed for better results (e.g. filtering bike lanes, stairs, etc.).

In order to be used for real routing / isochrone computation, we need to:

• improve snapping to edges (closest point, not vertex)
• split the edges at intersections (pgr_separateTouching / pgr_separateCrossing)
• handle turn restrictions
• split the road into smaller segments to have better speed representation

All can be achieved with pgRouting/PostGIS functions, a good dataset and time.

links

https://docs.pgrouting.org/3.8/en
https://qgis.org/