# Changes between Version 9 and Version 10 of UsersWikiSimplifyPreserveTopology

Show
Ignore:
Timestamp:
04/08/12 14:45:16 (14 months ago)
Comment:

--

Unmodified
Removed
Modified
• ## UsersWikiSimplifyPreserveTopology

v9 v10
164164}}}
165165
166== Function wrapper ==
167
168The 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):
169
170{{{
171-- takes a table name, gid and geom column and simplify the layer based on given tolerance.
172-- returns a setof (gid int, geom geometry). where geom is the simplified geometry and gid is the tablename gid.
173create or replace function simplifyLayerPreserveTopology (tablename text, id_col text, geom_col text, tolerance float) returns setof record as \$\$
174    DECLARE
175        tabname alias for \$1;
176        tid alias for \$2;
177        geo alias for \$3;
178        tol alias for \$4;
179        numpoints int:=0;
180
181    BEGIN
182
183        EXECUTE 'select sum(st_npoints(geom)) from '||quote_ident(tabname) into numpoints;
184        raise notice 'Num points in %: %', tabname, numpoints;
185
186        EXECUTE 'create unlogged table poly as ('
187                ||'select '||quote_ident(tid)||', (st_dump('||quote_ident(geo)||')).* from '||quote_ident(tabname)||')';
188
189        -- extract rings out of polygons
190        create unlogged table rings as
191        select st_exteriorRing((st_dumpRings(geom)).geom) as g from poly;
192
193        raise notice 'rings created...';
194
195        drop table poly;
196
197        -- 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
202        drop table rings;
203
204        -- extract lines as individual objects, in order to rebuild polygons from these
205        -- simplified lines
206        create unlogged table simplelines as select (st_dump(g)).geom as g from simplerings;
207
208        drop table simplerings;
209
210        -- Rebuild the polygons, first by polygonizing the lines, with a
211        -- distinct clause to eliminate overlaping segments that may prevent polygon to be created,
212        -- then dump the collection of polygons into individual parts, in order to rebuild our layer.
213        create temp table simplepolys as
214        select (st_dump(st_polygonize(distinct g))).geom as g
215        from simplelines;
216
217        select count(*) from simplepolys into numpoints;
218
219        raise notice 'rings polygonized. num rings: %', numpoints;
220
221        drop table simplelines;
222
223        -- some spatial indexes
224        create index simplepolys_geom_gist on simplepolys  using gist(g);
225
226        raise notice 'spatial index created...';
227
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)
230        RETURN QUERY EXECUTE 'select '||quote_ident(tid)||', st_collect(geom) as geom '
231            ||'from ('
232            --||'    select distinct on (d.'||quote_ident(tid)||') d.'||quote_ident(tid)||', s.g as geom '
233            ||'    select d.'||quote_ident(tid)||', s.g as geom '
234            ||'   from '||quote_ident(tabname)||' d, simplepolys s '
235            ||'    where st_intersects(d.'||quote_ident(geo)||', s.g) '
236            ||'    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 '
238            ||'    ) as foo '
239            ||'group by '||quote_ident(tid);
240
241        drop table simplepolys;
242
243        RETURN;
244
245    END;
246\$\$ language plpgsql strict;
247}}}
248
249=== Example usage ===
250
251{{{
252select * from simplifyLayerPreserveTopology('departement', 'gid', 'geom', 10000) as (gid int, geom geometry);
253}}}
254
255== ToDo ==
256
257Handle enclosed polygons cases.
166258Nicolas Ribot.