Changes between Version 9 and Version 10 of UsersWikiSimplifyPreserveTopology

Timestamp:

04/08/12 14:45:16 (14 months ago)

Author:

nicolasribotosgeo

Comment:

--

UsersWikiSimplifyPreserveTopology

@@ -164 +164 @@
 }}}
 
+== 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.