Splitting a table of linestring by a table of points onto or close to the lines, using Linear Referencing functions

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 lines | 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.

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 an ascending id for each successive point location.

Then, a self join of the table will be used to select the 2 locations to give to st_line_substring

 
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;

Splitting is very fast compare to other method using st_split(geom, st_expand(point))

Nicolas Ribot.