Changes between Version 1 and Version 2 of UsersWikiSimplifyPreserveTopology

Timestamp:

04/08/12 06:21:14 (13 years ago)

Author:

nicolasribotosgeo

Comment:

--

UsersWikiSimplifyPreserveTopology

@@ -3 +3 @@
 == What we want ==
 
-Simplifying this layer [[Image(SPT_dept_ori.png)]] into this: [[image:SPT_dept_sim.png]] 
+Simplifying this layer:
+
+ [[Image(SPT_dept_ori.png)]] 
+
+into this: 
+
+ [[Image(SPT_dept_sim.png)]] 
 
 == The data ==
@@ -24 +30 @@
 == Steps ==
 
+''Steps are divided into individual queries for the sake of clarity. The process could be done in a PL/PGSQL function to generalize it'''
+
 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.
 
 {{{
-create table dep_poly as select gid, code, st_dump
+create table poly as (
+    select gid, code_dept, (st_dump(geom)).* 
+    from departement
+);
 }}} 
+
+2. extract rings out of polygons
+
+{{{
+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:
+
+{{{
+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
+
+{{{
+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. 
+{{{
+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
+
+{{{
+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.
+
+{{{
+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. It gives better results in our testcase.
+{{{
+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
+);
+}}}
+
+9. rebuild departements by grouping them by code_dept (other attributes could be re-associated here): 
+
+{{{
+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)]]