# Changes between Version 16 and Version 17 of UsersWikiSimplifyPreserveTopology

Show
Ignore:
Timestamp:
04/09/12 11:03:58 (14 months ago)
Comment:

--

Unmodified
Removed
Modified
• ## UsersWikiSimplifyPreserveTopology

v16 v17
105105
106 8. Second method is based on percentage of overlaping area comparison.
107 {{{
108 #!sql
1068. Second method is based on percentage of overlaping area comparison. Empirical ratio used here.
107{{{
108#!sql
109
109110create table simpledep as (
110         select distinct on (d.code_dept) d.code_dept, s.g as geom
111        select d.code_dept, s.g as geom
111112        from departement d, simplepolys s
112113        where st_intersects(d.geom, s.g)
113         order by d.code_dept, st_area(st_intersection(s.g, d.geom))/st_area(s.g) desc
114        and st_area(st_intersection(s.g, d.geom))/st_area(s.g) > 0.5
114115);
115116}}}

148149        select gid, code_dept, (st_dump(geom)).*
149150        from departement
150 ) select distinct on (d.code_dept) d.code_dept, baz.geom
151) select d.code_dept, baz.geom
151152 from (
152153        select (st_dump(st_polygonize(distinct geom))).geom as geom

161162poly d
162163where st_intersects(d.geom, baz.geom)
163 order by d.code_dept, st_area(st_intersection(baz.geom, d.geom))/st_area(baz.geom) desc;
164and st_area(st_intersection(d.g, baz.geom))/st_area(baz.g) > 0.5;
164165}}}
165166

178179        tol alias for \$4;
179180        numpoints int:=0;
181        time text:='';
180182
181183    BEGIN
182184
183         EXECUTE 'select sum(st_npoints(geom)) from '||quote_ident(tabname) into numpoints;
184         raise notice 'Num points in %: %', tabname, numpoints;
185        EXECUTE 'select sum(st_npoints(geom)), to_char(clock_timestamp(), ''MI:ss:MS'') from '||quote_ident(tabname) into numpoints, time;
186        raise notice 'Num points in %: %. Time: %', tabname, numpoints, time;
185187
186188        EXECUTE 'create unlogged table poly as ('

191193        select st_exteriorRing((st_dumpRings(geom)).geom) as g from poly;
192194
193         raise notice 'rings created...';
195        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
196        raise notice 'rings created: %', time;
194197
195198        drop table poly;
196199
197200        -- Simplify the rings. Here, no points further than 10km:
198         create unlogged table simplerings as select st_simplifyPreserveTopology(st_linemerge(st_union(g)), tol) as g from rings;
199
200         raise notice 'rings simplified...';
201        create unlogged table gunion as select st_union(g) as g from rings;
202
203        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
204        raise notice 'union done: %', time;
201205
202206        drop table rings;
207
208        create unlogged table mergedrings as select st_linemerge(g) as g from gunion;
209
210        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
211        raise notice 'linemerge done: %', time;
212
213        drop table gunion;
214
215        create unlogged table simplerings as select st_simplifyPreserveTopology(g, tol) as g from mergedrings;
216
217
218        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
219        raise notice 'rings simplified: %', time;
220
221        drop table mergedrings;
203222
204223        -- extract lines as individual objects, in order to rebuild polygons from these

216235
217236        select count(*) from simplepolys into numpoints;
218
219         raise notice 'rings polygonized. num rings: %', numpoints;
237        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
238        raise notice 'rings polygonized. num rings: %. time: %', numpoints, time;
220239
221240        drop table simplelines;

226245        raise notice 'spatial index created...';
227246
228         -- comparing percentage of overlaping area gives good results to identify simplified objects.
229         -- (the 50% of overlapping area is empirical. Should fail with some odd-shaped polygons)
247        -- works better: comparing percentage of overlaping area gives better results.
248        -- as input set is multipolygon, we first explode multipolygons into their polygons, to
249        -- be able to find islands and set them the right departement code.
230250        RETURN QUERY EXECUTE 'select '||quote_ident(tid)||', st_collect(geom) as geom '
231251            ||'from ('

235255            ||'    where st_intersects(d.'||quote_ident(geo)||', s.g) '
236256            ||'    and st_area(st_intersection(s.g, d.'||quote_ident(geo)||'))/st_area(s.g) > 0.5 '
237             ||'    order by d.'||quote_ident(tid)||', st_area(st_intersection(s.g, d.'||quote_ident(geo)||'))/st_area(s.g) desc '
238257            ||'    ) as foo '
239258            ||'group by '||quote_ident(tid);