= An example showing how to simplify a multipolygon layer, keeping topology between objects =

== What we want ==

|| Simplifying this layer: || into this: || avoiding that: 
|| [[Image(SPT_dept_ori.png)]] ||  [[Image(SPT_dept_sim.png)]] || [[Image(SPT_no_topo.png)]] ||

== 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 ==

{{{
#!sh
shp2pgsql -IiD -g geom -s 2154 DEPARTEMENT.SHP departement | psql
}}}

== Principle of simplification ==

Based on the technique described here [http://gis.stackexchange.com/questions/178/simplifying-adjacent-polygons]:
 * 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''

1. 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.

{{{
#!sql
create table poly as (
    select gid, code_dept, (st_dump(geom)).* 
    from departement
);
}}} 

2. extract rings out of polygons

{{{
#!sql
create table rings as (
    select st_exteriorRing((st_dumpRings(geom)).geom) as g 
    from poly
);
}}}

3. 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:

{{{
#!sql
create table simplerings as (
    select st_simplifyPreserveTopology(st_linemerge(st_union(g)), 10000) as g 
    from rings
);
}}}

4. extract lines as individual objects, in order to rebuild polygons from these simplified lines

{{{
#!sql
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. 
{{{
#!sql
create table simplepolys as ( 
    select (st_dump(st_polygonize(distinct g))).geom as g
    from simplelines
);
}}}

6. Add an id column to help us identify objects and a spatial index

{{{
#!sql
alter table simplepolys  add column gid serial primary key;
create index simplepolys_geom_gist on simplepolys  using gist(g);
}}}

7. 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.

{{{
#!sql
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:

 [[Image(SPT_bad_attributes.png)]] 

8. Second method is based on percentage of overlaping area comparison.
{{{
#!sql
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:

 [[Image(SPT_good_attributes.png)]] 


9. rebuild departements by grouping them by code_dept (other attributes could be re-associated here): 

{{{
#!sql
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:

[[Image(SPT_simple_dept.png)]] 

Small island where collapsed during process:

[[Image(SPT_islands_removed.png)]] 

== One big query ==

{{{
#!sql
with poly as (
	select gid, code_dept, (st_dump(geom)).* 
	from departement
) select distinct on (d.code_dept) d.code_dept, baz.geom 
 from ( 
	select (st_dump(st_polygonize(distinct geom))).geom as geom
	from (
		select (st_dump(st_simplifyPreserveTopology(st_linemerge(st_union(geom)), 10000))).geom as geom
		from (
			select st_exteriorRing((st_dumpRings(geom)).geom) as geom
			from poly
		) as foo
	) as bar
) as baz,
poly d
where st_intersects(d.geom, baz.geom)
order by d.code_dept, st_area(st_intersection(baz.geom, d.geom))/st_area(baz.geom) desc;
}}}

== Function wrapper ==

The following code could be wrapped into a function like this. (It uses unlogged and temp tables, but not sure if it really helps to get speed):

{{{
-- takes a table name, gid and geom column and simplify the layer based on given tolerance.
-- returns a setof (gid int, geom geometry). where geom is the simplified geometry and gid is the tablename gid.
create or replace function simplifyLayerPreserveTopology (tablename text, id_col text, geom_col text, tolerance float) returns setof record as $$
    DECLARE
        tabname alias for $1;
        tid alias for $2;
        geo alias for $3;
        tol alias for $4;
        numpoints int:=0;

    BEGIN

        EXECUTE 'select sum(st_npoints(geom)) from '||quote_ident(tabname) into numpoints;
        raise notice 'Num points in %: %', tabname, numpoints;
        
        EXECUTE 'create unlogged table poly as ('
                ||'select '||quote_ident(tid)||', (st_dump('||quote_ident(geo)||')).* from '||quote_ident(tabname)||')';

        -- extract rings out of polygons
        create unlogged table rings as 
        select st_exteriorRing((st_dumpRings(geom)).geom) as g from poly;
        
        raise notice 'rings created...';
        
        drop table poly;

        -- Simplify the rings. Here, no points further than 10km:
        create unlogged table simplerings as select st_simplifyPreserveTopology(st_linemerge(st_union(g)), tol) as g from rings;
        
        raise notice 'rings simplified...';
        
        drop table rings;

        -- extract lines as individual objects, in order to rebuild polygons from these
        -- simplified lines
        create unlogged table simplelines as select (st_dump(g)).geom as g from simplerings;
        
        drop table simplerings;

        -- 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 temp table simplepolys as 
        select (st_dump(st_polygonize(distinct g))).geom as g
        from simplelines;
        
        select count(*) from simplepolys into numpoints;
        
        raise notice 'rings polygonized. num rings: %', numpoints;
        
        drop table simplelines;

        -- some spatial indexes
        create index simplepolys_geom_gist on simplepolys  using gist(g);

        raise notice 'spatial index created...';

        -- comparing percentage of overlaping area gives good results to identify simplified objects.
        -- (the 50% of overlapping area is empirical. Should fail with some odd-shaped polygons)
        RETURN QUERY EXECUTE 'select '||quote_ident(tid)||', st_collect(geom) as geom '
            ||'from ('
            --||'    select distinct on (d.'||quote_ident(tid)||') d.'||quote_ident(tid)||', s.g as geom '
            ||'    select d.'||quote_ident(tid)||', s.g as geom '
            ||'   from '||quote_ident(tabname)||' d, simplepolys s '
            ||'    where st_intersects(d.'||quote_ident(geo)||', s.g) '
            ||'    and st_area(st_intersection(s.g, d.'||quote_ident(geo)||'))/st_area(s.g) > 0.5 '
            ||'    order by d.'||quote_ident(tid)||', st_area(st_intersection(s.g, d.'||quote_ident(geo)||'))/st_area(s.g) desc '
            ||'    ) as foo '
            ||'group by '||quote_ident(tid);
            
        drop table simplepolys;
        
        RETURN;
    
    END;
$$ language plpgsql strict;
}}}

=== Example usage ===

{{{
select * from simplifyLayerPreserveTopology('departement', 'gid', 'geom', 10000) as (gid int, geom geometry);
}}} 

== !ToDo ==

Handle enclosed polygons cases.
Nicolas Ribot.