2

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: radars table

network table

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

3
  • Try breaking down that command into its components. For instance, just try SELECT ST_Line_Locate_Point(roads.geom, ST_ClosestPoint(roads.geom, radars.geom)) AS LR FROM... It's hard to tell what's wrong otherwise. It could be topology, it could be something in the code that I've overlooked. I just copied that code from an old file and changed it slightly so it had more general names. Which parts of the statement work, and which don't? Try it without each of the WHERE clauses, for example. Keep making it smaller. May 17 '15 at 9:58
  • When running that small statement (SELECT ST_Line_Locate_Point(network.geom, ST_ClosestPoint(network.geom, radars.geom)) AS LR FROM network, radars), I get the error - ERROR: line_locate_point: 1st arg isnt a line CONTEXT: SQL function "st_line_locate_point" statement 2 ********** Error ********** ERROR: line_locate_point: 1st arg isnt a line SQL state: XX000. So could the error be that it's a MultiLineString and not a LineString? May 17 '15 at 10:12
  • Yes, I think that is quite likely. You probably need singlepart lines for your network to make any sense. There's a tool to convert your lines that you can run first. May 17 '15 at 18:50
0

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);
0

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;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.