Version 6 (modified by 12 years ago) ( diff ) | ,
---|
An example showing how to simplify a multipolygon layer, keeping topology between objects
What we want
Simplifying this layer: | into this: |
![]() | ![]() |
The data
French administrative subdivisions, called "départements", will be used. Data can be downloaded here: http://http://professionnels.ign.fr/DISPLAY/000/528/175/5281750/GEOFLADept_FR_Corse_AV_L93.zip
Data projection is Lambert-93, EPSG:2154
Loading the data
shp2pgsql -IiD -g geom -s 2154 DEPARTEMENT.SHP departement | psql
Principle of simplification
Based on the technique described here [[]]:
- we extract linestrings out of polygons,
- then union and simplify them
- st_polygonize() is used to rebuild surfaces from linestrings.
- Finally, attributes from the initial layer are associated with simplified polygons.
Steps
Steps are divided into individual queries for the sake of clarity. One big query summarize them at the end of this page
- First extract the input multipolygons into polygons, keeping their departement code. This will allow us to associate attributes to each part of multipolygons at the end of the process.
create table poly as ( select gid, code_dept, (st_dump(geom)).* from departement );
- extract rings out of polygons
create table rings as ( select st_exteriorRing((st_dumpRings(geom)).geom) as g from poly );
- Simplify the rings. At this step, we choose the simplification ratio we want (some trials can be made by calling st_simplifyPreserveTopology on departement table).
Here, no points further than 10km:
create table simplerings as ( select st_simplifyPreserveTopology(st_linemerge(st_union(g)), 10000) as g from rings );
- extract lines as individual objects, in order to rebuild polygons from these simplified lines
create table simplelines as ( select (st_dump(g)).geom as g from simplerings );
5.rebuild the polygons, first by polygonizing the lines, with a distinct clause to eliminate overlaping segments that may prevent polygon to be created, then dump the collection of polygons into individual parts, in order to rebuild our layer.
create table simplepolys as ( select (st_dump(st_polygonize(distinct g))).geom as g from simplelines );
- Add an id column to help us identify objects and a spatial index
alter table simplepolys add column gid serial primary key; create index simplepolys_geom_gist on simplepolys using gist(g);
- attribute association between input layer and simplified polygons: First method to retrieve attribute is based on containment of a point of the surface of simplified polygons.
create table simpledep as ( select code_dept, g from departement d, simplepolys s where st_contains(d.geom, st_pointOnSurface(s.g)) );
It does not work: code_dept=92 is a curved polygon and fails with st_contains() test:
- Second method is based on percentage of overlaping area comparison.
create table simpledep as ( select distinct on (d.code_dept) d.code_dept, s.g as geom from departement d, simplepolys s where st_intersects(d.geom, s.g) order by d.code_dept, st_area(st_intersection(s.g, d.geom))/st_area(s.g) desc );
It gives better results in our testcase:
- rebuild departements by grouping them by code_dept (other attributes could be re-associated here):
create table simple_departement as ( select code_dept, st_collect(geom) as geom from simpledep group by code_dept );
Result
Layer now looks like this:
Small island where collapsed during process:
One big query
with poly as ( select gid, code_dept, (st_dump(geom)).* from departement ) select distinct on (d.code_dept) d.code_dept, qux.geom from ( select (st_dump(st_polygonize(distinct baz.geom))).geom as geom from ( select (st_dump(geom)).geom as geom from ( select st_simplifyPreserveTopology(st_linemerge(st_union(geom)), 10000) as geom from ( select st_exteriorRing((st_dumpRings(geom)).geom) as geom from poly ) as foo ) as bar ) as baz ) as qux, poly d where st_intersects(d.geom, qux.geom) order by d.code_dept, st_area(st_intersection(qux.geom, d.geom))/st_area(qux.geom) desc;
Nicolas Ribot.
Attachments (9)
-
SPT_dept_ori.png
(50.2 KB
) - added by 12 years ago.
Departement, originals
-
SPT_bad_attributes.png
(55.3 KB
) - added by 12 years ago.
simplified departements, bad attribute association
-
SPT_dept_sim.png
(45.6 KB
) - added by 12 years ago.
Final result
-
SPT_good_attributes.png
(55.2 KB
) - added by 12 years ago.
simplified departements, correct attribute association
-
SPT_islands_removed.png
(59.3 KB
) - added by 12 years ago.
simplified departements, zoom on simplified islands
-
SPT_simple_dept.png
(56.5 KB
) - added by 12 years ago.
simplified departements, with attributes
-
SPT_no_topo.png
(80.5 KB
) - added by 12 years ago.
simplified departements with broken topology
- world_before.png (72.0 KB ) - added by 12 years ago.
- world_after.png (67.4 KB ) - added by 12 years ago.
Download all attachments as: .zip