[GIS] Cutting linestrings with points

openstreetmappostgispostgresqlsplittingstreets

I was checking the best way to cut linestrings by points.

The scenario is: lots of streets, need the segments cut by intersect points, kind of this:

example

I got

  • linestrings (full uncuted by points) table

  • st_intersection points table

I need to get the independent linestring segments cut by the intersection points table.

I'm using PostGIS functions, and found several approaches, but every one of them has giving me some kind of issue.

Here's what I already tested:

1

Line table: 1 row , st_memunion of 1200 lines
Point table: 1700 rows (points)

What's bad: really takes really lot of time and memory flush. Can't create multiple tables at same time cause memory just can't handle it.
And the result is dirty and messy. Instead of giving me the correct row number, and I need to clean it up later (well explained here Split lines at intersection points)

CREATE TABLE lines_with_messy_result AS (
SELECT
((ST_DUMP(ST_SPLIT(a.geom,b.ix))).geom) as geom
FROM st_union_lines a
INNER JOIN lots_of_points b
ON ST_INTERSECTS(a.geom, b.ix)
);

--then need to clean this up
create table lines_segments_cleaned as (
SELECT DISTINCT ON (ST_AsBinary(geom)) 
geom 
FROM 
lines_with_messy_result
);

source of this way/approach: https://stackoverflow.com/questions/25753348/how-do-i-divide-city-streets-by-intersection-using-postgis

2

Same lines/points table.
Still messy results and need to clean this up.
Still lots of time to finish query.

--messy table    
Create table messy_lines_segments as 
Select 
row_number() over() as id,
(st_dump(st_split(input.geom, blade.ix))).geom as geom
from st_union_lines input
inner join lots_of_points blade on st_intersects(input.geom, blade.ix);

--then cleaning the messy table
delete from messy_lines_segments a
where exists 
(
select 1 from messy_lines_segments b where a.id != b.id and st_coveredby(b.geom,a.geom)
);

source of this way / approach : Split lines at intersection points

3

I also found this function:
https://github.com/Remi-C/PPPP_utilities/blob/master/postgis/rc_Split_Line_By_Points.sql

which has the good thing that it does not leave a messy_result that then I need to clean it up. But you do need st_memunion from both sides (lines table and points table)

It's kind of:

create table osm.calles_cortadas_03segmentos_sanluis as (
SELECT result.geom
FROM 
osm.calles_cortadas_01uniones_sanluis AS line, 
osm.calles_cortadas_00intersecciones_sanluis AS point,
rc_split_line_by_points(
input_line:=line.geom
,input_points:=point.ix
,tolerance:=4
) AS result
);

But its also super-long hours for getting the results. And I also tried with longer tables (10k lines, 14k points) and I just get memory out issues.

I also tried Esri's ArcGIS with also bad results…

So, what's the best way of getting this done with PostGIS geom functions?

I mean, without stepping into topology.

Or what's your best recommendation?

Best Answer

Best non-topology solution ever.. it REALLY works fast and easy (believe me i tested lots of ways to do this):

--please add this function:
https://github.com/Remi-C/PPPP_utilities/blob/master/postgis/rc_split_multi.sql

-- please create a universal sequence for unique id multipurpose
CREATE SEQUENCE select_id
  INCREMENT 1
  MINVALUE 1
  MAXVALUE 9223372036854775807
  START 1
  CACHE 1
  CYCLE;



-- lines (non st_unioned linestrings, with gist index) : linetable
-- points (non st_unioned points, with gist index) : pointtable

-- first, lines with cuts (pointtable)
create table segment_lines_00 as (
WITH points_intersecting_lines AS( 
   SELECT lines.id AS lines_id, lines.geom AS line_geom,  ST_Collect(points.ix) AS blade
   FROM 
   linetable as lines, 
   pointtable as points
   WHERE st_dwithin(lines.geom, points.ix, 4) = true
   GROUP BY lines.id, lines.geom
)
SELECT lines_id, rc_Split_multi(line_geom, blade, 4)
FROM points_intersecting_lines
);
CREATE INDEX ON segment_lines_00 USING gist(rc_split_multi);


-- then, segments cutted by points
create table segment_lines_01 as (
select 
(ST_Dump(rc_split_multi)).geom as geom,  nextval('SELECT_id'::regclass) AS gid  from segment_lines_00
);
CREATE INDEX ON segment_lines_01 USING gist(geom);
Related Question