Changes between Version 16 and Version 17 of UsersWikiSimplifyPreserveTopology

Timestamp:

Apr 9, 2012, 11:03:58 AM (12 years ago)

Author:

nicolasribotosgeo

Comment:

—

UsersWikiSimplifyPreserveTopology

@@ -104 +104 @@
  [[Image(SPT_bad_attributes.png)]] 
 
-8. Second method is based on percentage of overlaping area comparison.
-{{{
-#!sql
+8. Second method is based on percentage of overlaping area comparison. Empirical ratio used here.
+{{{
+#!sql
+
 create table simpledep as (
-        select distinct on (d.code_dept) d.code_dept, s.g as geom
+        select 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
+        and st_area(st_intersection(s.g, d.geom))/st_area(s.g) > 0.5
 );
 }}}
@@ -148 +149 @@
         select gid, code_dept, (st_dump(geom)).* 
         from departement
-) select distinct on (d.code_dept) d.code_dept, baz.geom 
+) select d.code_dept, baz.geom 
  from ( 
         select (st_dump(st_polygonize(distinct geom))).geom as geom
@@ -161 +162 @@
 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;
+and st_area(st_intersection(d.g, baz.geom))/st_area(baz.g) > 0.5;
 }}}
 
@@ -178 +179 @@
         tol alias for $4;
         numpoints int:=0;
+        time text:='';
 
     BEGIN
 
-        EXECUTE 'select sum(st_npoints(geom)) from '||quote_ident(tabname) into numpoints;
-        raise notice 'Num points in %: %', tabname, numpoints;
+        EXECUTE 'select sum(st_npoints(geom)), to_char(clock_timestamp(), ''MI:ss:MS'') from '||quote_ident(tabname) into numpoints, time;
+        raise notice 'Num points in %: %. Time: %', tabname, numpoints, time;
         
         EXECUTE 'create unlogged table poly as ('
@@ -191 +193 @@
         select st_exteriorRing((st_dumpRings(geom)).geom) as g from poly;
         
-        raise notice 'rings created...';
+        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
+        raise notice 'rings created: %', time;
         
         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...';
+        create unlogged table gunion as select st_union(g) as g from rings;
+        
+        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
+        raise notice 'union done: %', time;
         
         drop table rings;
+        
+        create unlogged table mergedrings as select st_linemerge(g) as g from gunion;
+        
+        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
+        raise notice 'linemerge done: %', time;
+        
+        drop table gunion;
+        
+        create unlogged table simplerings as select st_simplifyPreserveTopology(g, tol) as g from mergedrings;
+        
+        
+        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
+        raise notice 'rings simplified: %', time;
+        
+        drop table mergedrings;
 
         -- extract lines as individual objects, in order to rebuild polygons from these
@@ -216 +235 @@
         
         select count(*) from simplepolys into numpoints;
-        
-        raise notice 'rings polygonized. num rings: %', numpoints;
+        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
+        raise notice 'rings polygonized. num rings: %. time: %', numpoints, time;
         
         drop table simplelines;
@@ -226 +245 @@
         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)
+        -- works better: comparing percentage of overlaping area gives better results.
+        -- as input set is multipolygon, we first explode multipolygons into their polygons, to
+        -- be able to find islands and set them the right departement code.
         RETURN QUERY EXECUTE 'select '||quote_ident(tid)||', st_collect(geom) as geom '
             ||'from ('
@@ -235 +255 @@
             ||'    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);