Context Navigation

Changes between Version 2 and Version 3 of UsersWikiSplitPolygonWithPoints

Timestamp:: 04/27/13 14:13:49 (12 years ago)
Author:: nicolasribotosgeo
Comment:: --

Legend:

: Unmodified
: Added
: Removed
: Modified

UsersWikiSplitPolygonWithPoints

-              v2
+              v3
 Linear Referencing function st_line_interpolate_point() and st_line_substring() will be used to first locate all points along their respective Linestrings, then cut linestrings based on these locations.
+In the first CTE, an union is performed to generate the 0 and 1 line fractions, that are necessary to correctly split linestrings.
 To identify each point location for a Linestring, we use the rank() windowing function to generate a ascending id for each successive point location.
+Then, a self join of the table will be used to select 2 locations
+Based on the technique described here [http://gis.stackexchange.com/questions/178/simplifying-adjacent-polygons]:
+ * we extract linestrings out of polygons,
+ * then union and simplify them
+ * st_polygonize() is used to rebuild surfaces from linestrings.
+ * Finally, attributes from the initial layer are associated with simplified polygons.
+== Steps ==
+''Steps are divided into individual queries for the sake of clarity. One big query summarize them at the end of this page''
+. 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.
+Then, a self join of the table will be used to select the 2 locations to give to st_line_substring
 {{{
 #!sql
+create table poly as (
+    select gid, code_dept, (st_dump(geom)).*
+    from departement
+);
+}}}
+. extract rings out of polygons
+{{{
+#!sql
+create table rings as (
+    select st_exteriorRing((st_dumpRings(geom)).geom) as g
+    from poly
+);
+with locus as (
+        select l.gid, 0 as l
+        from lines l
+        UNION ALL
+        select l.gid, 1 as l
+        from lines l
+        UNION ALL
+        select l.gid, st_line_locate_point(l.geom, p.geom) as  l
+        from lines l, points p
+        where st_dwithin(l.geom, p.geom, 0.01)
+        order by gid, l
+), loc_with_idx as (
+        select gid, l, rank() over (partition by gid order by l) as idx
+        from locus
+)
+select l.gid, st_line_substring(l.geom, loc1.l, loc2.l) as geom
+from loc_with_idx loc1 join loc_with_idx loc2 using (gid) join lines l using (gid)
+where loc2.idx = loc1.idx+1;
 }}}
-. 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:
-{{{
-#!sql
-create table simplerings as (
-    select st_simplifyPreserveTopology(st_linemerge(st_union(g)), 10000) as g
-    from rings
-);
-}}}
-. extract lines as individual objects, in order to rebuild polygons from these simplified lines
-{{{
-#!sql
-create table simplelines as (
-    select (st_dump(g)).geom as g
-    from 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.
-{{{
-#!sql
-create table simplepolys as (
-    select (st_dump(st_polygonize(distinct g))).geom as g
-    from simplelines
-);
-}}}
-. Add an id column to help us identify objects and a spatial index
-{{{
-#!sql
-alter table simplepolys  add column gid serial primary key;
-create index simplepolys_geom_gist on simplepolys  using gist(g);
-}}}
-. 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.
-{{{
-#!sql
-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)]]
-. Second method is based on percentage of overlaping area comparison. Empirical ratio used here.
-{{{
-#!sql
-create table simpledep as (
-        select d.code_dept, s.g as geom
-        from departement d, simplepolys s
-        where st_intersects(d.geom, s.g)
-        and st_area(st_intersection(s.g, d.geom))/st_area(s.g) > 0.5
-);
-}}}
- It gives better results in our testcase:
- [[Image(SPT_good_attributes.png)]]
-. rebuild departements by grouping them by code_dept (other attributes could be re-associated here):
-{{{
-#!sql
-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)]]
-== One big query ==
-{{{
-#!sql
-with poly as (
-        select gid, code_dept, (st_dump(geom)).*
-        from departement
-) select d.code_dept, baz.geom
- from (
-        select (st_dump(st_polygonize(distinct geom))).geom as geom
-        from (
-                select (st_dump(st_simplifyPreserveTopology(st_linemerge(st_union(geom)), 10000))).geom as geom
-                from (
-                        select st_exteriorRing((st_dumpRings(geom)).geom) as geom
-                        from poly
-                ) as foo
-        ) as bar
-) as baz,
-poly d
-where st_intersects(d.geom, baz.geom)
-and st_area(st_intersection(d.g, baz.geom))/st_area(baz.g) > 0.5;
-}}}
-== 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):
-{{{
--- Simplify the given table of multipolygon with the given tolerance.
--- This function preserves the connection between polygons and try to avoid generating gaps between objects.
--- To identify objects after simplification, area comparison is performed, instead of PIP test, that may fail
--- with odd-shaped polygons. Area comparison may also failed on some cases
--- Example: (table 'departement' is in the public schema)
--- select * from simplifyLayerPreserveTopology('', 'departement', 'gid', 'geom', 10000) as (gid int, geom geometry);
---
--- @param schename: text, the schema name of the table to simplify. set to null or empty string to use search_path-defined schemas
--- @param tablename: text, the name of the table to simplify
--- @param idcol: text, the name of a unique table identifier column. This is the gid returned by the function
--- @param tolerance: float, the simplify tolerance, in object's unit
--- @return a setof (gid, geom) where gid is the identifier of the multipolygon, geom is the simplified geometry
-create or replace function simplifyLayerPreserveTopology (schemaname text, tablename text, idcol text, geom_col text, tolerance float)
-returns setof record as $$
-    DECLARE
-        schname alias for $1;
-        tabname alias for $2;
-        tid alias for $3;
-        geo alias for $4;
-        tol alias for $5;
-        numpoints int:=0;
-        time text:='';
-        fullname text := '';
-    BEGIN
-        IF schname IS NULL OR length(schname) = 0 THEN
-            fullname := quote_ident(tabname);
-        ELSE
-            fullname := quote_ident(schname)||'.'||quote_ident(tabname);
-        END IF;
-        raise notice 'fullname: %', fullname;
-        EXECUTE 'select sum(st_npoints('||quote_ident(geo)||')), to_char(clock_timestamp(), ''MI:ss:MS'') from '
-            ||fullname into numpoints, time;
-        raise notice 'Num points in %: %. Time: %', tabname, numpoints, time;
-        EXECUTE 'create unlogged table public.poly as ('
-                ||'select '||quote_ident(tid)||', (st_dump('||quote_ident(geo)||')).* from '||fullname||')';
-        -- extract rings out of polygons
-        create unlogged table rings as
-        select st_exteriorRing((st_dumpRings(geom)).geom) as g from public.poly;
-        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 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
-        -- 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.
-        drop table if exists simplepolys;
-        create  table simplepolys as
-        select (st_dump(st_polygonize(distinct g))).geom as g
-        from simplelines;
-        select count(*) from simplepolys into numpoints;
-        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
-        raise notice 'rings polygonized. num rings: %. time: %', numpoints, time;
-        drop table simplelines;
-        -- some spatial indexes
-        create index simplepolys_geom_gist on simplepolys  using gist(g);
-        raise notice 'spatial index created...';
-        -- 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('||quote_ident(geo)||') 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 '||fullname||' d, simplepolys s '
-            --||'    where (st_intersects(d.'||quote_ident(geo)||', s.g) or st_contains(d.'||quote_ident(geo)||', s.g))'
-            ||'    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 '
-            ||'    ) as foo '
-            ||'group by '||quote_ident(tid);
-        --drop table simplepolys;
-        RETURN;
-    END;
-$$ language plpgsql strict;
-}}}
-=== Example usage ===
-''(table 'departement' is in public schema)''
-{{{
-select * from simplifyLayerPreserveTopology('', 'departement', 'gid', 'geom', 10000) as (gid int, geom geometry);
-}}}
-== Example with countries ==
-Simplification tolerance: 0.3 degrees :)
-{{{
-create table simple_countries as (
-    select * from simplifyLayerPreserveTopology('', 'countries', 'gid', 'geom', 0.3) as (gid int, geom geometry);
-}}}
-|| [[Image(world_before.png)]] || [[Image(world_after.png)]] ||
-== !ToDo ==
-Handle enclosed polygons cases.