Topology cleaning with PostGIS

An early tester of the new PostGIS Topology submitted an interesting dataset which kept me busy for a couple of weeks fixing a bunch of bugs related to numerical stability/robustness.

Finally, the ST_CreateTopoGeo function succeeded and imported the dataset as a proper topological schema. Here’s what it looks like:

Edges of the built topology

At a first glance it doesn’t seem to be particularly problematic. Here’s the composition summary:

=# select topologysummary('small_sample_topo');
                    topologysummary
--------------------------------------------------------
 Topology small_sample_topo (2042), SRID 0, precision 0
 83 nodes, 156 edges, 74 faces, 0 topogeoms in 0 layers

But the devil hides at high zoom levels. Where to zoom ? What are we looking for ?

We are guaranteed none of the constructed edges cross so the only leftover problem we might encounter is very small faces constructed wherever the original input had small overlaps or underlaps (gaps). We can have a visual signal of those faces by creating a view showing faces with an area below a given threshold. Let’s do that:

CREATE VIEW small_sample_topo.small_areas AS
SELECT face_id, st_getfacegeometry('small_sample_topo', face_id)
FROM small_sample_topo.face
WHERE face_id > 0
AND st_area(st_getfacegeometry('small_sample_topo', face_id)) < 0.1;

That query would let us see where to find faces with area < 0.1 units. And here’s qgis showing it:

Areas smaller than 0.1 square units

Now we know where to zoom, and also the ID of the offending faces.

Let’s zoom in and show some labels:

Detail of small area

You can now see that face 59 is bound by (among others) edges 130 and 129. Just get rid of one of them to assign the area to an adjacent face.

We drop edge 130 using ST_RemEdgeModFace, assigning the area to face 52. Here’s the result:

Area after cleanup

Cleaning further would require removing further edges and thus getting rid of all the small faces. There’s a lot of room for automating such processes. The good new is you can now build your own automation around your specific use cases while still using robust and standard foundations.

It is to be noted that the whole process I described here only involved the geometrical/topological level and didn’t affect at all the semantic/feature level. If we had TopoGeometry objects defined by the faces we’d also know which small faces were part of overlaps or underlaps and could then act consequently by adding or removing faces to the definition of the appropriate TopoGeometry object. Such step would have been required for the overlap situations as the ST_RemEdgeModFace function doesn’t let you change the shape of a defined TopoGeometry.

Unfortunately the semantic level is lost when using the ISO functions, as the whole ISO topology model doesn’t deal with features at all. This is why I think PostGIS would benefit from having a function that converts your simple features into topologically-defined features by adding any missing primitive to a topology schema and constructing the feature for you. Such function, is only waiting for sponsorship to become a reality of PostGIS 2.0. If you like what we’re building for your data integrity, please consider supporting the effort!

Tags: , ,

6 Responses to “Topology cleaning with PostGIS”

  1. strk says:

    An addendum for the braves.

    Remove the longest edge bounding faces with area < 0.1: WITH small_faces AS ( SELECT face.face_id FROM small_sample_topo.face WHERE face.face_id > 0
    AND st_area(st_getfacegeometry(‘small_sample_topo’::character varying, face.face_id)) < 0.1 ), alldata AS ( SELECT f.face_id, e.edge_id, st_length(e.geom), max(st_length(e.geom)) over (partition by f.face_id) as mlen FROM small_sample_topo.edge_data e, small_faces f WHERE e.left_face = f.face_id or e.right_face = f.face_id ) SELECT st_remedgemodface('small_sample_topo', edge_id) FROM alldata WHERE st_length = mlen; Drop dangling edges: SELECT st_remedgemodface('small_sample_topo', edge_id) FROM small_sample_topo.edge_data WHERE left_face = right_face; Remove isolated nodes: SELECT st_remisonode('small_sample_topo', node_id) FROM small_sample_topo.node WHERE containing_face is not null;

  2. Chico says:

    First of all i gotta say great work i’m using Postgis Topologies for quite a lot of use cases.

    While looking through the functions code I was wondering if you could help me with a graph route network problem.

    The problem is I want to create a topology on a street network but only creating nodes and edges at noded intersections (not just a geometrical intersection so ST_Touches) to solve bridge and tunnel problems.

    Therefore do I need to edit the topogeo_addlinestring function where ST_Intersection is implemented?

    Thanks in advance

    Chico

  3. strk says:

    The PostGIS topology model follows the “Topology-Geometry” ISO model, which is not appropriate for route networks.
    Another ISO model “Topology-Network” is what you are after, which is not implemented in the core right now.
    Consider giving pg_routing a try. Or start work on a core implementation, it’d be fun ! 🙂

  4. Chico says:

    🙂 Yeah but pg_routing in particular doesn’t help here as it assumes network topologies are already exisiting (apart from import tools like osm2pgrouting et al. which is not satisfying for a few reasons) Alright will work my way through the code…

  5. Samuel says:

    Is there a function that would allow me to identify the nodes of a given face? I would like to programmatically detect sliver polygons by identifying which faces possess a small area, and then using ST_Azimuth to determine which of these contain very narrow vertices.

    Your process is brilliant, however the dataset I’m dealing with contains highly variable sizes (geology units) and shapes. The area of valid faces is thus sometimes much smaller than the area of a gigantic sliver, the only common denominator is the value of the angle between edges.

    Any insight on how to retrieve specific nodes?

  6. strk says:

    You can obtain the list of nodes bounding a specific face with a query like this:

    WITH edges AS (
    SELECT start_node, end_node
    FROM MYTOPO.edge
    WHERE left_face = MYFACEID
    OR right_face = MYFACEID
    )
    SELECT start_node FROM edges UNION
    SELECT end_node FROM edges;

    For what you need it would probably be faster to start the scan from the nodes, analyze angles of incident edges and report the faces of the offending edge pairs.

    Would make a good general function to add to PostGIS.

    See topology.getNodeEdges if you’re looking for inspiration.

Leave a Reply