Accueil / Blog / Métier / 2013 / Use PostGIS topologies to clean-up road networks

Use PostGIS topologies to clean-up road networks

Par Mathieu Leplatre publié 03/07/2013, édité le 21/02/2019
This article gives a few basics to get started with using the PostGIS topology extension

This article gives a few basics to get started with using the PostGIS topology extension.

We will take avandtage of topologies to clean-up a real topological road network, coming from Toulouse OpenData files.

Topological

A topology is a general concept, where objects are defined by their relationships instead of their geometries. Instead of lines, we manipulate edges, vertices and faces : might remind you the core concepts of graph theory.

A topological road network is supposed to have their lines (edges) connected at single points (nodes).

In this example dataset, JOSM validator detects not less than 1643 errors :) Broken connections, crossing lines ...

Let's clean this up !

Installation

On Ubuntu 12.04, you just have to install PostGIS :

sudo apt-add-repository -y ppa:ubuntugis/ppa
sudo apt-get update
sudo apt-get install -y postgresql postgis

The topology extension is installed by default. Just activate it in your database:

CREATE DATABASE "roadsdb";
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
SET search_path = topology,public;

Data Import

Load your shapefile (using command-line) like usual :

schema="public."
db="roadsdb"
user="postgres"
password="postgres"
host="localhost"
ogr2ogr -f "PostgreSQL" PG:"host=${host} user=${user} dbname=${db} password=${password}" -a_srs "EPSG:2154" -nln ${schema}roads -nlt MULTILINESTRING ROAD_SHAPEFILE.SHP

Create and associate the PostGIS topology:

SELECT topology.CreateTopology('roads_topo', 2154);
SELECT topology.AddTopoGeometryColumn('roads_topo', 'public', 'roads', 'topo_geom', 'LINESTRING');

Convert linestrings to vertices and edges within the topology :

-- Layer 1, with 1.0 meter tolerance
UPDATE roads SET topo_geom = topology.toTopoGeom(wkb_geometry, 'roads_topo', 1, 1.0);

From now on, we have a topology, whose imperfections were corrected. It smoothly merged all dirty junctions, whose defects were at most 1.0 meter wide.

You may enconter insertion problems : the tool fails [1] and aborts the whole transaction. Use this snippet to skip errors and go on with the next records:

DO $$DECLARE r record;
BEGIN
  FOR r IN SELECT * FROM roads LOOP
    BEGIN
      UPDATE roads SET topo_geom = topology.toTopoGeom(wkb_geometry, 'roads_topo', 1, 1.0)
      WHERE ogc_fid = r.ogc_fid;
    EXCEPTION
      WHEN OTHERS THEN
        RAISE WARNING 'Loading of record % failed: %', r.ogc_fid, SQLERRM;
    END;
  END LOOP;
END$$;

This is rather frustrating to face topological errors at insertion ! You can try with a lower tolerance, or check that your records have at least valid geometries. Any clarification or help on this would be welcome :)

Visualize and export

In order to visualize your topology vertices in QGIS, browse your database tables, and add the following layers: roads_topo.edge_data and roads_topo.node.

You can also export the resulting geometries into a new table :

CREATE TABLE roads_clean AS (
    SELECT ogc_fid, topo_geom::geometry
    FROM roads
);

Or obtain your lovable Shapefile in return :

ogr2ogr -f "ESRI Shapefile" ROAD_CLEAN.SHP PG:"host=${host} user=${user} dbname=${db} password=${password}" -sql "SELECT topo_geom::geometry FROM roads"

If, like Amit you want to split the lines at intersections and assign original attributes, just join roads_topo.edge_data and on the roads table :

SELECT r.lib_off, r.ogc_fid, e.geom
FROM roads_topo.edge e,
     roads_topo.relation rel,
     roads r
WHERE e.edge_id = rel.element_id
  AND rel.topogeo_id = (r.topo_geom).id

Going further...

We could collapse crossing lines and disconnected junctions into a nice and clean network.

Yes ahem, we weren't able to repair every topological error of this dataset using this automatic method. Some inconsistencies, like the following one, are like 6 meters wide ! They are, by the way, perfectly described in OpenStreetMap :

We could also play with simplifications using Sandro Santilli's SimplifyEdgeGeom [2] function, it will collapse edges with a higher tolerance ...

SELECT SimplifyEdgeGeom('roads_topo', edge_id, 1.0) FROM roads_topo.edge;

Don't hesitate to share your thoughts and feedback. Concrete use cases and examples are rare about this! And as usual, drop a comment if anything is wrong or not clear :)

[1] SQL/MM Spatial exception, geometry intersects edge, side location conflict, ...
[2] Just execute the function SQL code. It's just an elegant wrapper around ST_ChangeEdgeGeom and ST_Simplify.
ABONNEZ-VOUS À LA NEWSLETTER !
Voir aussi
Formation PostgreSQL / PostGIS du 26 au 28 mars à Paris Formation PostgreSQL / PostGIS du 26 au 28 mars à Paris 07/02/2019

Découvrez les outils Libres pour gérer vos données spatiales !

Détecter des formes dans des photos de paysage Détecter des formes dans des photos de paysage 14/02/2019

Et parce que c'est la St-Valentin, on détecte des cœurs !

Représentation des modèles numériques de terrain sur le web : ombrage et 3D Représentation des modèles numériques de terrain sur le web : ombrage et 3D 11/02/2019

Les Modèles Numériques de Terrain sont des données représentant la forme du terrain. Sur une ...

Cartographier le manteau neigeux avec Python Cartographier le manteau neigeux avec Python 08/01/2019

En quelques lignes de code, nous vous proposons de cartographier le manteau neigeux. Nous ...

Retour State of the Map France 2018 Retour State of the Map France 2018 23/07/2018

L'équipe de Makina attendait avec impatience l'événement de l'année, State Of The Map 2018 à ...