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:

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

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:

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.

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:

8. Second method is based on percentage of overlaping area comparison. Empirical ratio used here.

   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:

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:

Small island where collapsed during process:

One big query

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);

ToDo

Handle enclosed polygons cases.