Fork me on GitHub

Use PostGIS topologies to clean-up road networks

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.

Comments !

social