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