Tools: [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 ## 7. Visualization with points ## Conclusion 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) The result should be understood as a network reach envelope, not a real travel-time model. Upload road network to the database (gdal, qgis). The data used here wha downloaded from the french IGN https://geoservices.ign.fr/bdtopo ST_Force2D is used to drop Z/M values if present. Always add an index on the geometry column for performance. Add required columns for pgRouting. Replace NULL speeds with 0 (unusable edges). Compute costs based on length and speed (in km/h). With BD TOPO, I exaggerate the speed 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 Costs are expressed in seconds. Intersections are free (if a node exists ofc), speeds are approximate, results are optimistic. This is the point layer representing all start/end points of edges. Always add indexes for performance, well not always but when it makes sense. For each edge (line), find the start and end point and link to the corresponding vertex id (calculated previously). The result is visible in QGIS. cost and reverse_cost are in seconds. 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. 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. As we can see, the closest vertex is not very representative of the actual road network location. 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. And here is the complete query to compute the isochrones. Quick tip: As polygons are not donuts. You can still use QGIS to make sure the minutes taken are visible: 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. For each point in bdtopo.fire_stations, compute travel-time isochrones for: Using a road network stored in bdtopo.edges. Costs are expressed in seconds. Vertices are classified into exclusive time ranges: Each reachable node is joined to bdtopo.edges_vertices to retrieve its POINT geometry. For each (fire_station_id, minutes) group: These polygons represent: Approximate envelopes of reachable network vertices They are suitable for: They are not precise travel-time surfaces. 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. 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. You can visualize the result in QGIS. This helps you understand what can be wrong/explained with/by the network data. These isochrones are approximations. They are useful for exploration, not precision routing; unless you have a very good road network :) In order to be used for real routing / isochrone computation, we need to: All can be achieved with pgRouting/PostGIS functions, a good dataset and time. https://docs.pgrouting.org/3.8/en https://qgis.org/ Templates let you quickly answer FAQs or store snippets for re-use. Are you sure you want to hide this comment? It will become hidden in your post, but will still be visible via the comment's permalink. Hide child comments as well For further actions, you may consider blocking this person and/or reporting abuse COMMAND_BLOCK: 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 Enter fullscreen mode Exit fullscreen mode COMMAND_BLOCK: 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 COMMAND_BLOCK: 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 CODE_BLOCK: CREATE EXTENSION postgis; CREATE EXTENSION pgrouting; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: CREATE EXTENSION postgis; CREATE EXTENSION pgrouting; CODE_BLOCK: CREATE EXTENSION postgis; CREATE EXTENSION pgrouting; CODE_BLOCK: 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; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: 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; CODE_BLOCK: 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; CODE_BLOCK: CREATE INDEX edges_gix ON bdtopo.edges USING GIST (geom); Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: CREATE INDEX edges_gix ON bdtopo.edges USING GIST (geom); CODE_BLOCK: CREATE INDEX edges_gix ON bdtopo.edges USING GIST (geom); CODE_BLOCK: ALTER TABLE bdtopo.edges ADD COLUMN source bigint, ADD COLUMN target bigint, ADD COLUMN cost double precision, ADD COLUMN reverse_cost double precision; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: ALTER TABLE bdtopo.edges ADD COLUMN source bigint, ADD COLUMN target bigint, ADD COLUMN cost double precision, ADD COLUMN reverse_cost double precision; CODE_BLOCK: ALTER TABLE bdtopo.edges ADD COLUMN source bigint, ADD COLUMN target bigint, ADD COLUMN cost double precision, ADD COLUMN reverse_cost double precision; CODE_BLOCK: UPDATE bdtopo.edges SET vitesse_moyenne_vl = 0 WHERE vitesse_moyenne_vl IS NULL; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: UPDATE bdtopo.edges SET vitesse_moyenne_vl = 0 WHERE vitesse_moyenne_vl IS NULL; CODE_BLOCK: UPDATE bdtopo.edges SET vitesse_moyenne_vl = 0 WHERE vitesse_moyenne_vl IS NULL; CODE_BLOCK: 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); Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: 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); CODE_BLOCK: 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); CODE_BLOCK: CREATE TABLE bdtopo.edges_vertices AS SELECT * FROM pgr_extractVertices( $$SELECT id, geom FROM bdtopo.edges ORDER BY id$$ ); Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: CREATE TABLE bdtopo.edges_vertices AS SELECT * FROM pgr_extractVertices( $$SELECT id, geom FROM bdtopo.edges ORDER BY id$$ ); CODE_BLOCK: CREATE TABLE bdtopo.edges_vertices AS SELECT * FROM pgr_extractVertices( $$SELECT id, geom FROM bdtopo.edges ORDER BY id$$ ); CODE_BLOCK: CREATE INDEX edges_vertices_gix ON bdtopo.edges_vertices USING GIST (geom); CREATE INDEX edges_vertices_id ON bdtopo.edges_vertices (id); Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: CREATE INDEX edges_vertices_gix ON bdtopo.edges_vertices USING GIST (geom); CREATE INDEX edges_vertices_id ON bdtopo.edges_vertices (id); CODE_BLOCK: CREATE INDEX edges_vertices_gix ON bdtopo.edges_vertices USING GIST (geom); CREATE INDEX edges_vertices_id ON bdtopo.edges_vertices (id); CODE_BLOCK: 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; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: 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; CODE_BLOCK: 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; CODE_BLOCK: ALTER TABLE bdtopo.fire_stations ADD COLUMN IF NOT EXISTS start_vid bigint, ADD COLUMN IF NOT EXISTS snap_dist_m double precision; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: ALTER TABLE bdtopo.fire_stations ADD COLUMN IF NOT EXISTS start_vid bigint, ADD COLUMN IF NOT EXISTS snap_dist_m double precision; CODE_BLOCK: ALTER TABLE bdtopo.fire_stations ADD COLUMN IF NOT EXISTS start_vid bigint, ADD COLUMN IF NOT EXISTS snap_dist_m double precision; COMMAND_BLOCK: 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 ); Enter fullscreen mode Exit fullscreen mode COMMAND_BLOCK: 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 ); COMMAND_BLOCK: 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 ); CODE_BLOCK: CREATE TABLE bdtopo.isochrones ( fire_stations_id bigint, minutes int, geom geometry(Polygon, 2154) ); Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: CREATE TABLE bdtopo.isochrones ( fire_stations_id bigint, minutes int, geom geometry(Polygon, 2154) ); CODE_BLOCK: CREATE TABLE bdtopo.isochrones ( fire_stations_id bigint, minutes int, geom geometry(Polygon, 2154) ); CODE_BLOCK: 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; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: 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; CODE_BLOCK: 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; CODE_BLOCK: -- 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; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: -- 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; CODE_BLOCK: -- 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; CODE_BLOCK: (fire_station_id, node_id, travel_time_seconds) Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: (fire_station_id, node_id, travel_time_seconds) CODE_BLOCK: (fire_station_id, node_id, travel_time_seconds) CODE_BLOCK: 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; Enter fullscreen mode Exit fullscreen mode CODE_BLOCK: 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; CODE_BLOCK: 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; - speed information in BD TOPO is not always correct, and optimistic, work is being done to make it more accurate. - intersections are not modelled correctly (see pgr_separateTouching and/or use pgr_separateCrossing) https://docs.pgrouting.org/latest/en/pgRouting-concepts.html#id65 - https://docs.pgrouting.org/latest/en/pgRouting-concepts.html#id65 - snapping is approximate (closest vertex, not closest point on edge) - roads are sometimes included entirely even if only partially reachable - https://docs.pgrouting.org/latest/en/pgRouting-concepts.html#id65 - click on layer properties → Symbology → Categorized - check Control feature rendering order - pgr_drivingDistance is executed once per SIS point - Maximum distance: 20 minutes (1200 seconds) - Output nodes are filtered to agg_cost <= 20*60 - Points are collected - A tight concave hull is computed using alpha 0.10 - A small buffer (5 m) smooths the geometry - exploratory analysis - accessibility screening - visual communication - 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