Use PostGIS topologies to clean-up road networks
source link: https://blog.mathieu-leplatre.info/use-postgis-topologies-to-clean-up-road-networks.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
Use PostGIS topologies to clean-up road networks
Wed 03 July 2013This 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.#postgis, #sql, #topology, #opendata, #toulouse - Posted in the Dev category
© Copyright 2020 by Mathieu Leplatre. mnmlist Theme
Content licensed under the Creative Commons attribution-noncommercial-sharealike License.
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK