[[TOC]]

= 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:  
|| [[Image(SLBP_dept_ori.png)]] ||  [[Image(SLBP_points.png.png)]] 
||'''We want to split the linestrings by the points''' || 
|| [[Image(SLBP_res.png)]] ||

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

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

{{{
#!sql
 
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 compared with other method using st_split(geom, st_expand(point))

Nicolas Ribot.