Changes between Version 1 and Version 2 of UsersWikiSplitPolygonWithPoints

Timestamp:

Apr 27, 2013, 2:04:47 PM (11 years ago)

Author:

nicolasribotosgeo

Comment:

—

UsersWikiSplitPolygonWithPoints

@@ -1 @@
--- Splitting a table of linestring by a table of points onto or close to the lines --
+[[TOC]]
+
+= Splitting a table of linestring by a table of points onto or close to the lines  =
+
+== What we want ==
+
+|| Given this layer of Linestrings: || and this layer of points:  
+|| [[Image(SLBP_dept_ori.png)]] ||  [[Image(SLBP_points.png)]] 
+
+We want to split the linestrings by the points
+
+== The data ==
+
+French administrative subdivisions, called "départements", will be used. Data can be downloaded here: [http://professionnels.ign.fr/DISPLAY/000/528/175/5281750/GEOFLADept_FR_Corse_AV_L93.zip]
+The LIMITE_DEPARTEMENT.SHP shapefile is used. It contains departements limits as MultiLinestrings
+
+Data projection is Lambert-93, EPSG:2154
+
+== Loading the data ==
+
+{{{
+#!sh
+shp2pgsql -IiDS -g geom -s 2154 LIMITE_DEPARTEMENT.SHP departement | psql
+}}}
+
+The table of points will be generated by randomly choosing points on the linestrings:
+
+{{{
+#!sql
+create sequence points_seq;
+create table points as with rand as (
+        select random() as locus from generate_series (1, 10)
+), rand2 as (
+        select nextval('points_seq') as gid, l.gid as lgid, r.locus, st_line_interpolate_point(l.geom, locus) as geom
+        from rand as r, lines l
+) select * from rand2
+order by random() limit 1500; -- only 1500 pts out of 3300 generated with random locus on the linestrings
+
+}}}
+
+== Principle of Splitting ==
+
+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.
+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''
+
+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.
+
+{{{
+#!sql
+create table poly as (
+    select gid, code_dept, (st_dump(geom)).* 
+    from departement
+);
+}}} 
+
+2. extract rings out of polygons
+
+{{{
+#!sql
+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:
+
+{{{
+#!sql
+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
+
+{{{
+#!sql
+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. 
+{{{
+#!sql
+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
+
+{{{
+#!sql
+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.
+
+{{{
+#!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)]] 
+
+8. 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)]] 
+
+
+9. 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.