DEV Community: jwaddle

[PostGIS/PgRouting] With french routable data

jwaddle — Tue, 27 Jan 2026 12:10:14 +0000

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:

Points are collected
A tight concave hull is computed using alpha 0.10
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/

GTFS data powered by OpenTripPlanner - Part 1

jwaddle — Mon, 28 Feb 2022 09:55:38 +0000

GTFS data powered by OpenTripPlanner - Part 1

Introduction to GTFS and OTP.

What is GTFS?

The General Transit Feed Specification (GTFS) has been around for quite some time. First released 15 years ago GTFS started out as a side project of Google employees.
It can be used for network topology management, studying for a new network, new lines, or new timetables. But also for visualizing the topology of several networks, accessibility, analysis and having real-time information.

A GTFS feed is a collection of at least six, and up to 13 files contained within a zip file. These static text data files are consistent in format, incorporating timetable, geographic data and general agency data into a tabulated format (feed). GTFS “feeds” let public transit agencies publish the transit data.

Table name	Description	Sample
Agency	gathers transport service information (transport companies, network name)	agency_id: SEM agency_name: Mobilitées M -Tag agency_url: https://www.mobilites-m.fr/,Europe/Paris agency_timezone:null agency_lang:FR agency_phone:0438703870
Routes	Transit routes. A route is a group of trips that are displayed to riders as a single service.	agency_id: SEM route_id:1 route_short_name:C1 route_long_name:GRENOBLE Cité Jean Macé / MEYLAN Maupertuis / MONTBONNOT-SAINT- MARTIN Pré de l'Eau route_type:3 route_color:FDEA00 route_text_color: 000000
Trips	Trips for each route. A trip is a sequence of two or more stops that occur during a specific time period.
Stop times	Times that a vehicle arrives at and departs from stops for each trip.	trip_id; 25485431 stop_id; 6806 arrival_time; 07:22:00,07:22:00 departure_time; 07:22:00 Stop_sequence; 1 Pickup_type; 0
Calendar dates	Service dates specified using a weekly schedule with start and end dates. This file is required unless all dates of service are defined in calendar_dates.txt. And so, Exceptions for the services defined in the calendar.txt. If calendar.txt is omitted, then calendar_dates.txt is required and must contain all dates of service.	service_id: 12345-MFHLV13c_ST0_13_HI2122 date: 20220117 exception_type : 1
Transfer	Shows the correspondences between several breakpoints	from_stop_id: S6 to_stop_id: S7 transfer_type: 2 min_transfer_time: 300

Many transit agencies have created and published GTFS data with the primary purpose being integration with Google Maps. Data can be found fairly easily, the following website lists available feeds: https://transitfeeds.com/feeds

During testing we will be using data from the Grenoble area which can be found at data.metropolegrenoble.fr

What is OTP (Open Trip Planner)?

OTP (https://www.opentripplanner.org/) is an OpenSource Software written in Java, that is designed to use open data sources OpenStreetMap (OSM) and the GTFS data. It allows users to plan a trip that can combine multiple modes of transportation, such as bicycling or walking to reach public transportation. It can deal with Interlines (transfers) between different public transportation systems. Users can specify departure or arrival times for public transportation to build isochrones.

OTP is used in production in several cities, here a few examples

Alencon (France) : https://altobus.com/mon-itineraire
Atlanta (USA) : http://itsmarta.com/planatrip.aspx
Grenoble (France) : https://www.mobilites-m.fr/

OTP is available in two major versions that are still maintained. While the second version is optimized, it does not have all the functionalities from the first version, as demonstrated in the table below

TLDR

Open Trip Planner 1

Exists since 2009
Widespread use
Many options and tools
Transit routing approach obsolete
Experimental code

Pros of Open Trip Planner 2

2018
Still in Testing
Large or dense network
EU-standard Netex and SIRI data
More Efficient

Check http://docs.opentripplanner.org/en/latest/Version-Comparison/ for all differences

With no isochrones and no transfer policy with OTP 2 we opted
to use OTP 1 for the upcoming tutorial.

Use Cases
The following map (build with QGIS) shows our area of study with the list of stops from the GTFS data.

OTP provides several APIs (http://dev.opentripplanner.org/apidoc/1.5.0/index.html). We will cover a few of them.

The main feature of OpenTripPlanner is for planning trips. You can use the API to request an
itinerary from one place to another at a specific time:

http://0.0.0.0:8089/otp/routers/default/plan?fromPlace=45.16385418613824%2C5.732674598693848&toPlace=45.170117418455064%2C5.747051239013672&time=9%3A48am&date=01-28-2022&mode=TRANSIT%2CWALK&maxWalkDistance=750&arriveBy=false&wheelchair=false&locale=fr

key	value	description
fromPlace	45.16,5.73	the origin of the trip, in latitude, longitude
toPlace	45.17,5.74	the destination of the trip
time	9:48am	the desired departure time
date	01-28-2022	the desired departure date
mode	TRANSIT,WALK	transport modes to consider, in this case a combination of walking and transit
maxWalkDistance	750	the maximum distance in meters that you are willing to walk
arriveBy	false	specifies that the given time is when we plan to depart rather than when we want to arrive
wheelchair	false	If there is wheelchair access or not

This will give you a list of possible itineraries (as shown with the builtin website):

OpenTripPlanner can also be used to calculate the area which is accessible from a point within a given travel time, also known as a travel time isochrone.

GET http://localhost:8089/otp/routers/default/isochrone\?fromPlace\=45.18831,5.713\&mode\=WALK,TRANSIT\&maxWalkDistance\=500\&\&cutoffSec\=1000\&cutoffSec\=1800\&cutoffSec\=3600\&cutoffSec\=4800

HTTP 200
{
"result": "ok",
{
"type": "FeatureCollection",
"name": "export_4",
"crs": { "type": "name", "properties": {... } },
"features": [...], 
}

The cutoffSec param (1000, 3600…) tells the router that we want an isochrone ending at 1000, 1800, 3600, 4800 seconds) travel from the origin.

The previous request will result in a json file that once added to QGIS will show you the areas you can go to in 1000, 1800(seconds), ect.. .

Another powerful feature is the one-to-many travel time analysis. As an example, we want to find out how much time is needed to go to each available job using public transport.

Using the pointSet (a csv file with x, y representing jobs) and the surface API we are able to calculate the time required to go from one point (your home for example) to a potential job.
http://localhost:8089/otp/surfaces?batch=true&fromPlace=45.18831,5.713&time=5:46pm&date=01-20-2022&mode=TRANSIT,WALK&maxWalkDistance=800&arriveBy=false

Then querying http://localhost:8089/otp/surfaces/0/indicator?detail=true&targets=jobs&detail=true. Will return the amount of time needed to get to each job. Thus enabling us to generate the following map. Where each point represents a job and the label the number of minutes to get there.

That's it for this first part!
You should now be familiar with both GTFS data and the OTP application.

The next blog will continue with Part 2 and consist of a short tutorial in order to get you started with OTP.

Get in contact with us!
If you already want to receive additional information please do not hesitate to contact us!

DEV Community: jwaddle

[PostGIS/PgRouting] With french routable data

Road-network isochrones with PostGIS & pgRouting (Docker)

1. PostgreSQL + pgRouting in Docker

2. Preparing the Road Network

3. Travel Time (Length × Speed)

4. Graph Vertices

5. Snapping Origin Points

6. Isochrone Computation

Some explanation about the above query

SQL Workflow

Step-by-Step Explanation

1. Compute reachable vertices within 20 minutes (dd20)

2. Assign vertices to time bands (bucketed)

3. Convert vertices to points (reach_pts)

4. Build polygons from points

Interpretation

Tip

7. Visualization with points

Conclusion

links

GTFS data powered by OpenTripPlanner - Part 1

GTFS data powered by OpenTripPlanner - Part 1

Introduction to GTFS and OTP.

What is GTFS?

What is OTP (Open Trip Planner)?

1. Compute reachable vertices within 20 minutes (`dd20`)

2. Assign vertices to time bands (`bucketed`)

3. Convert vertices to points (`reach_pts`)