Splitting lines table using point table in PostGIS

Question

My problem is that I want to get the shortest route between multiple point pairs which are stored in a table (PostGIS).

In PostGIS I also have a network created using pgr_createTopology on a table containing about a 1000 MultiLineString rows. It's my understanding that to calculate routes between points with pg_routing, you need to either accept that the route goes to and from the nearest network node or you need to work around that limitation by adding new nodes to the network which corresponds to the stops you want to calculate on.

To add new nodes to my network I've tried tried to follow the answer from @alpha-beta-soup in this tread
- Split lines into non-overlapping subsets based on points
Before following the answer I used ArcMap to snap the points to the lines. I then imported the two shapefiles (network and radars both with geometry columns "geom") into PostGIS using the PostGIS Shapefile loader (2.1).

When using alpha_beta_soup's code I get a table with all the right columns, but there are zero rows in it, I've tried running only the first part of the code, but I still only get a table with no content. This is the first part of the code:

CREATE TABLE tmp_lineswithstops AS (
WITH subq AS (
    SELECT
    ST_Line_Locate_Point(
        roads.geom,
        ST_ClosestPoint(roads.geom, radars.geom)
    ) AS LR,
    rank() OVER (
        PARTITION BY roads.gid
        ORDER BY ST_Line_Locate_Point(
            roads.geom,
            ST_ClosestPoint(roads.geom, radars.geom)
        )
    ) AS LRRank,
    ST_ClosestPoint(roads.geom, radars.geom),
    roads.*
    FROM network AS roads
    LEFT OUTER JOIN radars
    ON ST_Distance(roads.geom, radars.geom) < 0.0001
    WHERE ST_Equals(ST_StartPoint(roads.geom), radars.geom) IS false
    AND ST_Equals(ST_EndPoint(roads.geom), radars.geom) IS false
    ORDER BY gid, LRRank
)
SELECT LR FROM subq
);  

The two tables look like this:

Does anyone have a clue as to what i'm doing wrong?

Answer 1

In case anyone lands here searching how to split lines with points, here is a current (2019) solution.

• union point geometry
• snap lines to points with tolerance of 1m - for projected coordinate system (if using geographic you'll need to adjust this value to degrees (0.00001 very roughly equals 1m)
• split lines at points
• insert source lines that weren't processed because they didn't intersect any points

id column must be unique or this won't work correctly:

create table edges_split as
  select a.id, (st_dump(st_split(st_snap(a.geom,b.geom,1),b.geom))).geom as geom from
    edges a, 
    (select id, st_union(geom) as geom from points) b;
insert into edges_split (id, geom) 
  select id, geom from edges 
  where id not in (select id from edges_split);

Answer 2

I had the same issue a number of years ago.

What worked for me is highlighted in part here: https://gis.stackexchange.com/a/272265/10661

I created a function for it, full code as a function below.

There was an issue splitting with points for me, what I found that worked consistently was converting the points to small lines that traverse the network and then splitting with those.

My methodology was as follows:

• Find nearest line to the input point(s) [called destination points below]
• Get a small subsection of this line
• Rotate the subsection 90 degrees around a point on the line closest to the input point (allows us to cut and or create point at intersection)
• Split the network table with the new lines [I called them blades]

/*
mc_cushroute_createsplitnetwork
A FUNCTION TO SPLIT A ROAD NETWORK AT DESTINATION POINTS AND CREATE A ROUTABLE NETOWRK FROM THAT.  
INPUT VARIBALES:  
    - SCHEMA NAME (name of the schema where the data sits - intermediate data will be created here also)  
    - NETWORK NAME (name of the network to split)  
    - DESTINATION NAME (name of the table containing destination points - these will be used to split the network [add nodes at destinatin locations])  
Michael Cushen    
@cushenmichael    
04/12/2017  
##########  
*/
--DROP FUNCTION mc_cushroute_createsplitnetwork(text,text,text);
CREATE OR REPLACE FUNCTION public.mc_cushroute_createsplitnetwork(schemaname text, networkname text, destinationname text)

    RETURNS TEXT AS

    $func$

    DECLARE

        a TEXT;

        b TEXT;

    BEGIN

        BEGIN

            a = 'DROP TABLE IF EXISTS '|| schemaname||'.test_closestline';

            EXECUTE a;
    a = 'CREATE TABLE '|| schemaname||'.test_closestline AS
        SELECT pid ,lid, distance_m, a.geom 
        FROM
            (SELECT q.id pid, a.id lid, ROUND(ST_Distance(q.geom, a.geom)::NUMERIC, 2) as distance_m, 
        ROW_NUMBER() OVER (PARTITION BY q.id ORDER BY ST_DISTANCE(a.geom,q.geom)) as row_number,
        a.geom
            FROM '|| schemaname|| '.' || networkname || ' a '
            'INNER JOIN '|| schemaname|| '.' ||  destinationname  || ' q ON
            ST_DWITHIN(a.geom,q.geom,100) -- use a sensible distance to restrict the result set without losing records 
            ) a
        WHERE a.row_number = 1';
    EXECUTE a;
RAISE NOTICE 'Closest lines created // test_closestline';   

    a = 'DROP TABLE IF EXISTS '|| schemaname||'.test_blade_rotate;';
    EXECUTE a;
    a = 'CREATE TABLE '|| schemaname||'.test_blade_rotate AS   
    SELECT
    line.id, ST_Rotate(line.geom, 1.5708, ST_Centroid(line.geom)) geom --1.5708 radians = 90 degrees
    FROM
       (
       SELECT paired.pid id, 
    CASE 
        WHEN  ( round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2) + 1 = 1 ) THEN  ST_LineSubstring(paired.lgeom, round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2), round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2)+ 0.01)
        WHEN  ( round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2) + 1 = 2) THEN  ST_LineSubstring(paired.lgeom, round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2) - 0.01, round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2))
        ELSE ST_LineSubstring(paired.lgeom,(round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2) - 0.01), (round((ST_LineLocatePoint(paired.lgeom,paired.pgeom))::numeric,2)+0.01))
    END AS geom 
       FROM 
            (SELECT line.pid, (ST_DUMP(line.geom)).geom lgeom, (ST_DUMP(point.geom)).geom pgeom
            FROM '
            || schemaname||'.test_closestline'||' as line,'
            || schemaname||'.' ||destinationname||' as point
            WHERE line.pid = point.id) as paired
       ) as line;';
    EXECUTE a;     
RAISE NOTICE 'Cutting blades created // test_blade_rotate'; 

    a = 'DROP INDEX IF EXISTS '|| schemaname||'.' || networkname || '_gix;';
    EXECUTE a;
    a = 'CREATE INDEX ' || networkname || '_gix ON '|| schemaname||'.' || networkname || ' USING GIST (geom);';
    EXECUTE a;
    a = 'DROP INDEX IF EXISTS '|| schemaname||'.test_blade_rotate_gix;';
    EXECUTE a;
    a = 'CREATE INDEX test_blade_rotate_gix ON '|| schemaname||'.test_blade_rotate USING GIST (geom);';
    EXECUTE a;
RAISE NOTICE 'indexes created'; 

    a = 'DROP TABLE IF EXISTS '|| schemaname||'.test_blade_split;';
    EXECUTE a;
    a = 'CREATE TABLE '|| schemaname||'.test_blade_split AS   
        SELECT line.id id, (ST_DUMP(line.geom)).geom  geom
        FROM 
        '|| schemaname||'.' || networkname || ' line
        WHERE   NOT EXISTS
            (
            SELECT  geom
            FROM    '|| schemaname||'.test_blade_rotate blade
            WHERE ST_INTERSECTS(line.geom, blade.geom)
            );'

    'DROP TABLE IF EXISTS '|| schemaname||'.test_blade_rotate_collect;'
    'CREATE TABLE '|| schemaname||'.test_blade_rotate_collect AS
    SELECT 1 id, ST_Collect(geom) geom
    FROM '|| schemaname||'.test_blade_rotate;'  

    'INSERT INTO '|| schemaname||'.test_blade_split
    SELECT cut.* 
    FROM
        (
        SELECT 
        line.id id, 
        (ST_DUMP(ST_SPLIT(line.geom, blade.geom))).geom geom
        FROM 
        '|| schemaname||'.' || networkname || ' line,
        '|| schemaname||'.test_blade_rotate_collect blade
        WHERE ST_INTERSECTS(line.geom, blade.geom)
        ORDER by line.id
        ) AS cut;';
    EXECUTE a;      
RAISE NOTICE 'Split Network created';   

--MIGHT WANT TO SPLIT NETWORK [ST_SPLIT] AT EACH JUNCTION AT THIS STAGE ALSO (WAS A PROBLEM WITH ARMADALE NETWORK)

    a = 'DROP SEQUENCE IF EXISTS '|| schemaname||'.test_blade_split_seq;';
    EXECUTE a;  
    a = 'CREATE SEQUENCE '|| schemaname||'.test_blade_split_seq
      START WITH 1
      INCREMENT BY 1
      MINVALUE 1
      NO MAXVALUE
      CACHE 1;';
    EXECUTE a;
    a = 'UPDATE '|| schemaname||'.test_blade_split set id = nextval('''|| schemaname||'.test_blade_split_seq'');';
    EXECUTE a;    
    a = 'DROP INDEX IF EXISTS '|| schemaname||'.test_blade_split_gix ;';
    EXECUTE a;  
    a = 'CREATE INDEX test_blade_split_gix ON '|| schemaname||'.test_blade_split USING GIST (geom);';
    EXECUTE a;  
RAISE NOTICE 'Unique ids added to Split Network';   


    a = 'ALTER TABLE '|| schemaname||'.test_blade_split ADD COLUMN "source" integer;';
    EXECUTE a;
    a = 'ALTER TABLE '|| schemaname||'.test_blade_split ADD COLUMN "target" integer;';
    EXECUTE a;
    a = 'DROP TABLE IF EXISTS '|| schemaname||'.test_blade_split_vertices_pgr;';
    EXECUTE a;
END;


RETURN 'ALL DONE :D';
END
$func$ LANGUAGE plpgsql;
--run function and create topology for split network
BEGIN;
SELECT mc_cushroute_createsplitnetwork('schema','routing_network','destination_points');
COMMIT;
ALTER TABLE schema.routing_network ALTER COLUMN id TYPE integer;
BEGIN;
SELECT pgr_createtopology('schema.test_blade_split', 0.00001, 'geom', 'id');
COMMIT;
BEGIN;
SELECT pgr_analyzeGraph('schema.test_blade_split', 0.000001, the_geom := 'geom', id := 'id');
COMMIT;
BEGIN;
ALTER TABLE schema.test_blade_split ADD COLUMN "cost" float8;
UPDATE schema.test_blade_split SET cost=ST_length(geom);
ALTER TABLE schema.test_blade_split_vertices_pgr  RENAME COLUMN the_geom TO geom;
COMMIT;