## Simplifying a map layer using PostGIS topology

Following a recent research about how to simplify a multipolygon layer while keeping topological relationships intact, here’s my take on that, using the PostGIS topological support.

### The data

It is composed by 96 multipolygons for a total of  47036 vertices.

### Principle of simplification

• We convert a layer’s Geometries to TopoGeometries
• We simplify all edges of the built topology
• We convert the (now-simplified) TopoGeometries back to Geometries

### Steps

The following steps assume you loaded the shapefile into a table named “france_dept”.

``` -- Create a topology
SELECT CreateTopology('france_dept_topo', find_srid('public', 'france_dept', 'geom'));

SELECT AddTopoGeometryColumn('france_dept_topo', 'public', 'france_dept', 'topogeom', 'MULTIPOLYGON');

-- Populate the layer and the topology
UPDATE france_dept SET topogeom = toTopoGeom(geom, 'france_dept_topo', 1); -- 8.75 seconds

-- Simplify all edges up to 10000 units
SELECT SimplifyEdgeGeom('france_dept_topo', edge_id, 10000) FROM france_dept_topo.edge; -- 3.86 seconds

-- Convert the TopoGeometries to Geometries for visualization
ALTER TABLE france_dept ADD geomsimp GEOMETRY;
UPDATE france_dept SET geomsimp = topogeom::geometry; -- 0.11 seconds```

### The SimplifyEdgeGeom function

You may have noticed that the “SimplifyEdgeGeom” is not a core function. It is a function I wrote for the purpose of catching topological problems introduced by simplification.

The naive call would be:

```SELECT ST_ChangeEdgeGeom('france_dept_topo', edge_id, ST_Simplify(geom, 10000))
FROM france_dept_topo.edge;```

The problem with the above call is that any simplification introducing a topology error would be rejected by ST_ChangeEdgeGeom by throwing an exception and the exception would rollback the whole transaction leaving you with no edge changed. Possible topology errors introduced are: edges collapsing to points, intersecting self or other edges.

The SimplifyEdgeGeom function wraps the ST_ChangeEdgeGeom call into a subtransaction and handles exceptions by reducing the simplification factor until it succeeds. The version I used reduces simplification factor in half at each failure, dropping down to zero around 1e-8. You can roll your own with other heuristics or generalize this one to take parameters about stepping and limits.

Here’s the function:

```CREATE OR REPLACE FUNCTION SimplifyEdgeGeom(atopo varchar, anedge int, maxtolerance float8)
RETURNS float8 AS \$\$
DECLARE
tol float8;
sql varchar;
BEGIN
tol := maxtolerance;
LOOP
sql := 'SELECT topology.ST_ChangeEdgeGeom(' || quote_literal(atopo) || ', ' || anedge
|| ', ST_Simplify(geom, ' || tol || ')) FROM '
|| quote_ident(atopo) || '.edge WHERE edge_id = ' || anedge;
BEGIN
RAISE DEBUG 'Running %', sql;
EXECUTE sql;
EXCEPTION
WHEN OTHERS THEN
RAISE WARNING 'Simplification of edge % with tolerance % failed: %', anedge, tol, SQLERRM;
tol := round( (tol/2.0) * 1e8 ) / 1e8; -- round to get to zero quicker
IF tol = 0 THEN RAISE EXCEPTION '%', SQLERRM; END IF;
END;
END LOOP;
END
\$\$ LANGUAGE 'plpgsql' STABLE STRICT;```

### Performance

The times shown near the “expensive” steps give you an indication of the performance you may expect. It’s about 13 seconds in total for the 3 steps outlined in the first paragraph.

A single run of the simplification step brought vertices down to 1369 (from 47036).

Timings are take on this system:

POSTGIS=”2.0.1SVN r9637″ GEOS=”3.4.0dev-CAPI-1.8.0″ PROJ=”Rel. 4.8.0, 6 March 2012″ GDAL=”GDAL 1.9.0, released 2011/12/29″ LIBXML=”2.7.6″ LIBJSON=”UNKNOWN” TOPOLOGY RASTER

PostgreSQL 8.4.10 on x86_64-pc-linux-gnu, compiled by GCC gcc-4.4.real (Ubuntu 4.4.3-4ubuntu5) 4.4.3, 64-bit

CUSTOM OPTIONS:
shared_buffers = 128MB      (default is 24MB)
temp_buffers = 32MB         (default is 8MB)
work_mem = 8MB              (default is 1MB)
maintenance_work_mem = 32MB (default is 16MB)
max_stack_depth = 4MB       (default is 2MB)
checkpoint_segments = 24    (default is 3)

CPU: Intel(R) Core(TM)2 Duo CPU P9500 @ 2.53GHz (2 x 5054.34 bogomips)

RAM: 4GB

### Considerations

The procedure described in this post is also valid for LINESTRING and MULTILINESTRING layers, using the exactly same code.

You could reuse the topology to produce multiple resolution levels w/out incurring again in the construction cost (and with a reduced input complexity at each level).

The simplification step doesn’t use TopoGeometry objects at all so you could choose to perform  topology construction and  attribute assignment in a different way.

Running the SimplifyEdgeGeom function again might give you more simplification because edges which may have intersected to the simplified version of an edge may not be intersecting anymore after their own simplification. The function can be changed to behave differently on exception to improve performance or quality.

Tags: ,

### 27 Responses to “Simplifying a map layer using PostGIS topology”

1. […] Simplifying a map layer using PostGIS topology « Strk's Blog. Publicado: 2012/04/13 Arquivado em: geo, postgis […]

2. James Crone says:

awesome post. Got round to installing postgis 2 this morning and was looking to try out generalising polygons using topology. This recipe really helped.

3. strk says:

An update: this exercise exposed a bug in ST_ChangeEdgeGeom accepting a theoretically invalid move. I figured because the move back was correctly identified as invalid (non isomorphic).

See more here: http://trac.osgeo.org/postgis/ticket/1775

The fix will be released in 2.0.1

4. Jacolin says:

Hello,

GeoFla files are now available with this URL: http://professionnels.ign.fr/sites/default/files/GEOFLADept_FR_Corse_AV_L93.zip (EPSG code is 2154).

Regards,

Y.

5. Phil Stephens (peas) says:

When I run the first command, I get the following error:

```ERROR: function createtopology(unknown, integer) does not exist LINE 1: SELECT CreateTopology('usqctsmod_topo', find_srid('public', ... ^ HINT: No function matches the given name and argument types. You might need to add explicit type casts.```

``` ********** Error ********** ```

```ERROR: function createtopology(unknown, integer) does not exist SQL state: 42883 Hint: No function matches the given name and argument types. You might need to add explicit type casts. Character: 8 ```

6. strk says:

Make sure you enabled topology: `select postgis_full_version();` should contain a TOPOLOGY label when enabled.
Also make sure your search_path contains the “topology” schema (topolog.CreateTopology would work if that’s the problem).

7. peas says:

select postgis_full_version(); returns this:
“POSTGIS=”2.0.0 r9605″ GEOS=”3.3.4-CAPI-1.7.3″ PROJ=”Rel. 4.8.0, 6 March 2012″ GDAL=”GDAL 1.9.1, released 2012/05/15″ LIBXML=”2.7.8″ LIBJSON=”UNKNOWN” TOPOLOGY RASTER”

I’m unaware as to how to check the search_path in any database. ?

8. peas says:

quick update: I was able to create a topology by calling topology.CreateTopology()

9. strk says:

Use: `SHOW search_path;` to get the current value and `SET search_path TO x,y,z` for setting the value.

10. Francis Seong says:

I want to reduce edges and nodes. I run “select st_remedgemodface(‘topo’, edge_id)”, I got follow error:

ERROR: TopoGeom #number in layer 1 (#table.topogeom) cannot be represented dropping edge #edge_id

11. strk says:

The message is telling you that by dropping that edge it will be impossible for one TopoGeometry object to retain the same shape.

You can see the shape of the TopoGeometry using:

`SELECT topogeom::geometry FROM #table WHERE id(topogeom) = #number`

12. Francis Seong says:

I run :
SELECT st_astext(topogeom::geometry) FROM idns_link_0 WHERE id(topogeom) = #table.id

result is :
“MULTILINESTRING((14128213.2930175 4505878.66922969,14128214.0420175 4505877.92022969))”

I don’t understant this result.

13. strk says:

I suggest subscribing to the postgis users mailing list for this kind of help

14. […] saw how to get a topologically consistent simplified version of a full layer, but that method isn’t fast enough for on-demand usage. Also, we saw how to perform a fast […]

15. Fabian Zeindl says:

when trying to use my own data, I get the following error:

UPDATE rayons SET topogeom = toTopoGeom(geom, ‘france_dept_topo’, 1);
ERROR: lwpoly_from_lwlines: shell must be closed
CONTEXT: SQL statement “WITH ids as ( select row_number() over () as seq, edge from unnest(‘{-160,-104,45,-46,41,-44,43,-53,50,-49,47,-67,66,64,-57,55,-61,60,-84,70,-73,72,79,-76,74,-102,101,-87,85,-96,95,-112,98,-100,93,-92,90,-121,120,-125,124,-119,117,-128,113,-109,107,115,-132,-131,-115,116,-108,114,-127,129,-39,-158,159,-153,155,-137,-138,147,-141,144,-145,152,-148,149,-133,134}’::int[] ) u(edge) ), edges AS ( select CASE WHEN i.edge < 0 THEN ST_Reverse(e.geom) ELSE e.geom END as g FROM ids i left join france_dept_topo.edge_data e ON(e.edge_id = abs(i.edge)) ORDER BY seq) SELECT ST_MakePolygon(ST_MakeLine(g.g)) FROM edges g;"
PL/pgSQL function _st_addfacesplit(character varying,integer,integer,boolean) line 34 at EXECUTE statement
SQL statement "SELECT topology._ST_AddFaceSplit(atopology, -newedge.edge_id, newedge.left_face, false)"
PL/pgSQL function st_addedgemodface(character varying,integer,integer,geometry) line 513 at SQL statement
PL/pgSQL function topogeo_addlinestring(character varying,geometry,double precision) line 124 at assignment
SQL statement "SELECT array_cat(edges, array_agg(x)) FROM ( select topology.TopoGeo_addLinestring(atopology, rec.geom, tol) as x ) as foo"
PL/pgSQL function topogeo_addpolygon(character varying,geometry,double precision) line 24 at assignment
PL/pgSQL function totopogeom(geometry,character varying,integer,double precision) line 100 at FOR over SELECT rows

any ideas?

16. strk says:

Please use the bug tracker for this: http://trac.osgeo.org/postgis

17. Wayne says:

I ran into an issue with an exception being thrown using the vanilla code in the blog post. I found that it was resolved by changing the order of the lines that calculate tol and the if condition that throws an error if tol is 0. This forces the code to try tolerance 0 first, then thrown an error if that fails.

SQL/MM Spatial exception is the exception that gets thrown.

I changed the code to the following:
IF tol = 0 THEN
RAISE EXCEPTION ‘%’, SQLERRM;
END IF;
tol := round( (tol/2.0) * 1e8 ) / 1e8; — round to get to zero quicker

Needless to say simplifying with tolerance 0 is not supported.

18. Michael says:

I’m trying to run this but I get the following error
psql:topology.sql:34: ERROR: relation “topology” does not exist

Line 34 is
UPDATE my_table SET topogeom = topology.toTopoGeom(geom, ‘my_table_topo’, 1);

The other functions work by prepending “topology.” but this doesn’t seem to do the trick here.

19. strk says:

Does the topology.topology table exist ? It should. Try \set VERBOSITY verbose

20. Benbannock says:

Super Post ! Merci

21. […] simplifier un ensemble de polygones limitrophes (départements) […]

22. If you receive an error like the following:

ERROR: function createtopology(unknown, integer) does not exist

you probably haven’t created the topology extension. You’ll have to run:

create extension postgis_topology;

(I also had to restart my psql prompt for the functions to register)

23. MatPas says:

Hi,
thank you for this feature.

I’m in a trouble trying to adapt topology to my needs, basically I need to delete portion of a layer to substitute it with a new part, this because I work on tiles.

I guess working with tiles should rose a lot of questions, but at this moment I’ll be very happy to find a way to delete a selection of edges and put there a square tile to stay in topology.

I’ve tried by building a list of edge_id starting from edge_data.geomXcut.aoi (or even passing throughout faces.geomXcut.aoi), retrieving which nodes and which faces are not implied in relations outside the geometry selection(cut.aoi), got the tree list (egdes, nodes, faces) the delete fail even inside the transaction (edge_data, const:initially deferrable property).

Can someone give an advise on this?

24. strk says:

You should split edges on the tile grid boundaries, which can be done by passing each grid line to TopoGeo_addLinestring.
The TopoGeometry definitions will be maintained to still cover the same pointset as they did before the cutting.

25. Gabriel Daldegan says:

Hi,
I’m trying to use topology to fill small and big gaps between polygons of very different sizes, and I’m getting the SQL/MM Spatial exception curve not simple SQL state message every time I try to run it with a tolerance higher than 4 units (meters). This tolerance is too small and is not filling the gaps I want.
Can anyone give a tip on how to proceed?

Thanks

26. strk says:

To fill the gaps you’ll want to first import the topology and then remove the faces which are smaller than a given tolerance.
Removing a face can be done by removing one of its bounding edges.

27. Gabriel Daldegan says:

Thank you Sandro!
I’ll try it doing as you suggested.