Implementing Point-based Costs in pgRouting

Implementing Point-based Costs in pgRouting

pgRouting is a great tool but it does have certain limitations.

One limitation is that one cannot have costs at specific nodes.

All costs have to be edge based.

Discussion on Node/Edge Costs in pgModeler

For example, suppose a network has a stop valve between two pipes. Normally the stop value is open (no cost to traverse) but when it is closed pgRouting should not be able to create a network trace that crosses the closed value.

Since pgRouting does not have the ability to create costs via points, and needs all costs to be expressed on edges, a method is needed to enable such a cost to be added to the network in order to affect any network tracing pgRouting may execute.

Here is my method for doing this.

(All scripts are available via a link at the end of the article.)

1. Create and Load a table containing network edges described by linestrings.

-- Create and load table.
--
DROP TABLE IF EXISTS public.network CASCADE;
CREATE TABLE public.network (
    id integer NOT NULL,
    material character varying(5),
    diameter double precision,
    zone_name double precision,
    pipeclass character varying(20),
    measured_length double precision,
    geom public.geometry(MultiLineString,28356)
);
	
INSERT INTO public.network (geom , id , material, diameter, zone_name, pipeclass) VALUES ('0105000020C46E000001000000010200000002000000BB490C824B461E41AE47E1CA2C815A4177BE9F1A4C461E41E17A149E2C815A41', 1, 'DI', 300, 29, 'UNK');
... etc ...

-- Set costs based on line length
--
UPDATE public.network SET measured_length = ST_Length(geom);

-- Create spatial index
--
CREATE INDEX network_geom_geom_idx ON public.network USING GIST (geom);

-- Add networking columns

ALTER TABLE public.network 
  ADD COLUMN source       integer,
  ADD COLUMN target       integer,
  ADD COLUMN cost         double precision DEFAULT 1,
  ADD COLUMN reverse_cost double precision DEFAULT 1;

vacuum analyze verbose public.network;
-- 3563 live rows

2. Create Stop Valves.

Now that the network table is created, a set of stop valves can be created.

The stop valves are created where a pipe crosses from one zone (region) to another. Zones group the pipes into management areas.

Where two connected pipes have different non-null zones, a valve will be created at the point at which they touch.

-- Create table to hold stop valves.
--
DROP TABLE IF EXISTS public.network_stops;

CREATE TABLE public.network_stops (
  id         serial,
  zones      varchar(100),
  node_point geometry(Point,28356)
);

-- Find connected pipes with different zone_names and insert them into the network_stops table.
--
INSERT INTO public.network_stops (zones,node_point) 
SELECT REPLACE(REPLACE(zones,',(null)',''),'(null),',''), node_point
  FROM (SELECT node_point,
               string_agg(distinct COALESCE(CAST(f.zone_name as varchar(100)),'(null)'),',') as zones
          FROM (SELECT zone_name, ST_SnapToGrid(ST_StartPoint(ST_GeometryN(geom,1)),0.001) as node_point 
                  FROM public.network AS a 
                UNION ALL
                SELECT zone_name, ST_SnapToGrid(ST_EndPoint(ST_GeometryN(geom,1)),0.001) as node_point 
                  FROM public.network AS b
               ) as f
         GROUP By node_point
       ) as g
 WHERE zones <> '(null)' 
   AND zones LIKE '%,%';

ALTER TABLE public.network_stops ADD CONSTRAINT network_stops_pk PRIMARY KEY (id);

CREATE INDEX network_stops_node_sx ON public.network_stops USING GIST (node_point);

ALTER TABLE public.network_stops CLUSTER ON network_stops_node_sx;

vacuum analyze verbose public.network_stops;
-- 28 live rows 

3. Add new short stop lines

At each stop add the short lines that fill the gaps left when the original input lines are shortened.

WITH data AS (
  -- Find all network lines entering/leaving a network stop
  SELECT wp.id,
         ST_GeometryN(wp.geom,1) as line,
         ns.node_point as point
    FROM public.network_stops as ns
         INNER JOIN public.network as wp
            ON ST_DWithin(wp.geom,ns.node_point,0.01)
)
, ShortLines as (
-- Shorten existing lines, and create short lines to fill gap.
SELECT id, 
       shortened_line, 
       short_line
  FROM (SELECT id,
               shortened_line, 
               CASE WHEN startEndType = 'START'
                    THEN ST_MakeLine(point,ST_StartPoint(shortened_line))
                    ELSE ST_MakeLine(ST_EndPoint(shortened_line),point)
                END as short_line
          FROM (SELECT id,
                       point, 
                       CASE WHEN ST_Equals(ST_StartPoint(line),point) 
                            THEN 'START'
                            ELSE 'END'
                        END as startEndType,
                       ST_GeometryN(
                            CASE WHEN ST_Equals(ST_StartPoint(line),point) 
                                 THEN spdba.ST_Reduce(line,0.05,'START',3,2) 
                                 ELSE spdba.ST_Reduce(line,0.05,'END',3,2)
                             END,
                            1
                       ) as shortened_line
                  FROM data as a
               ) as f
       ) as g
 WHERE ROUND(ST_LENGTH(short_line)::numeric,3) < 0.055
)

-- Insert short_lines representing network stops along with original attributes.
-- Note, initially, their costs are set to -1 ("thou shalt not pass")
--
INSERT INTO public.network (
  source, target, cost, reverse_cost, material, diameter, pipeclass, zone_name, measured_length, geom
)
SELECT nw.source, nw.target, 
       -1 as cost,
       -1 as reverse_cost, 
       nw.material, nw.diameter, nw.pipeclass, nw.zone_name,
       ROUND(ST_Length(a.short_line)::numeric,3) as measured_length
       ST_Multi(a.short_line)
  FROM ShortLines as a
        inner join 
          public.network as nw
            on (nw.id = a.id);

4. Shorten original lines

Finally, update the network to shorten the original lines entering/leaving network stop point.

WITH data AS (
  -- Find all network lines entering/leaving a network stop
  SELECT wp.id,
         ST_GeometryN(wp.geom,1) as line,
         ns.node_point as point
    FROM public.network_stops as ns
         INNER JOIN public.network as wp
            ON ST_DWithin(wp.geom,ns.node_point,0.01)
)
, ShortLines as (
-- Shorten existing lines, and create short lines to fill gap.
SELECT id, 
       shortened_line, 
       short_line
  FROM (SELECT id,
               shortened_line, 
               CASE WHEN startEndType = 'START'
                    THEN ST_MakeLine(point,ST_StartPoint(shortened_line))
                    ELSE ST_MakeLine(ST_EndPoint(shortened_line),point)
                END as short_line
          FROM (SELECT id,
                       point, 
                       CASE WHEN ST_Equals(ST_StartPoint(line),point) 
                            THEN 'START'
                            ELSE 'END'
                        END as startEndType,
                       ST_GeometryN(
                            CASE WHEN ST_Equals(ST_StartPoint(line),point) 
                                 THEN spdba.ST_Reduce(line,0.05,'START',3,2) 
                                 ELSE spdba.ST_Reduce(line,0.05,'END',3,2)
                             END,
                            1
                       ) as shortened_line
                  FROM data as a
               ) as f
       ) as g
 WHERE ROUND(ST_LENGTH(short_line)::numeric,3) < 0.055
)
-- Now update existing lines with their shortened equivalent
UPDATE public.network 
   SET geom = ST_Multi(ShortLines.shortened_line) 
  FROM (SELECT distinct id, shortened_line FROM ShortLines as s) as ShortLines
 WHERE network.id = ShortLines.id;
-- 56 Lines Updated

An example of a network stop and its new linestrings can be seen here.

Stop Lines
Network Stop with linear representation

5. Routing Examples.

5.1. Create routing network

-- First update public.topology to ensure it is fit for use for pgRouting
--
select pgr_createTopology('public.network',0.001,'geom','id','source','target',rows_where:='true',clean:=false);

/*
NOTICE:  PROCESSING:
NOTICE:  pgr_createTopology('public.network', 0.001, 'geom', 'id', 'source', 'target', rows_where := 'true', clean := f)
NOTICE:  Performing checks, please wait .....
NOTICE:  Creating Topology, Please wait...
NOTICE:  1000 edges processed
NOTICE:  2000 edges processed
NOTICE:  3000 edges processed
NOTICE:  -------------> TOPOLOGY CREATED FOR  3619 edges
NOTICE:  Rows with NULL geometry or NULL id: 0
NOTICE:  Vertices table for table public.network is: public.network_vertices_pgr
NOTICE:  ----------------------------------------------

Successfully run. Total query runtime: 5 secs 432 msec.
1 rows affected.
*/

-- Conduct necessary analysis
--
select pgr_analyzegraph('public.network', 0.001,the_geom:='geom',id:='id',source:='source',target:='target',rows_where:='true');

/*
NOTICE:  PROCESSING:
NOTICE:  pgr_analyzeGraph('public.network',0.001,'geom','id','source','target','true')
NOTICE:  Performing checks, please wait ...
NOTICE:  Analyzing for dead ends. Please wait...
NOTICE:  Analyzing for gaps. Please wait...
NOTICE:  Analyzing for isolated edges. Please wait...
NOTICE:  Analyzing for ring geometries. Please wait...
NOTICE:  Analyzing for intersections. Please wait...
NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 2
NOTICE:                            Dead ends: 342
NOTICE:  Potential gaps found near dead ends: 0
NOTICE:               Intersections detected: 24
NOTICE:                      Ring geometries: 0

Successfully run. Total query runtime: 1 secs 519 msec.
1 rows affected.
*/

5.2 Conduct a Trace With All Stops Open

Firstly we will execute a network trace without any restriction at the network stops

The cost/reverse_code for the ordinary lines was set when creating public.network.

-- Set "open" cost for stops
--
update public.network set cost = 1.1, reverse_cost = 1.1 where cost = -1;

-- Now run test between two points (see public.network_vertices_pgr) 
--
SELECT ST_AsEWKT(ST_Union(wp.geom)) as geom
  FROM pgr_dijkstra(edges_sql:='SELECT id, source, target, cost, reverse_cost FROM public.network',
                    start_vid:=(select id FROM public.network_vertices_pgr as a WHERE ST_DWithin(a.the_geom,'SRID=28356;POINT(498341.9 6948542.7)'::geometry,0.1)),
                      end_vid:=(select id FROM public.network_vertices_pgr as a WHERE ST_DWithin(a.the_geom,'SRID=28356;POINT(498560.5 6948744.7)'::geometry,0.1)),
                     directed:=false
       ) as f 
       INNER JOIN public.network as wp
          ON (wp.id = f.edge);

Visualising the result produces the following:

Open Stops
Open Stop 27

5.3 Conduct a Trace With Stop 27 Closed

We now set Stop to closed and then re-run the network trace.

-- Identify stop, find its in/out short lines, then update their costs
--
WITH data AS (
  SELECT nw.id,
         ns.node_point as point
    FROM public.network_stops as ns
           INNER JOIN public.network as nw
              ON ST_DWithin(nw.geom,ns.node_point,0.01)
  WHERE ns.id = (select id FROM public.network_vertices_pgr as a WHERE ST_DWithin(a.the_geom,'SRID=28356;POINT(498341.9 6948542.7)'::geometry,0.1))
    AND nw.cost = 1.1
)
UPDATE public.network 
   SET cost = -1, 
       reverse_cost = -1
  FROM data as d
 WHERE network.id = d.id
   AND cost = 1.1;

-- Re-run test
--
SELECT ST_AsEWKT(ST_Union(wp.geom)) as geom
  FROM pgr_dijkstra(edges_sql:='SELECT id, source, target, cost, reverse_cost FROM public.network',
                    start_vid:=(select id FROM public.network_vertices_pgr as a WHERE ST_DWithin(a.the_geom,'SRID=28356;POINT(498341.9 6948542.7)'::geometry,0.1)),
                      end_vid:=(select id FROM public.network_vertices_pgr as a WHERE ST_DWithin(a.the_geom,'SRID=28356;POINT(498560.5 6948744.7)'::geometry,0.1)),
                     directed:=false
       ) as f 
       INNER JOIN public.network as wp
          ON (wp.id = f.edge);

The result looks like this.

Trace after closing Stop 27
Closed Stop 27 trace.

6. Conclusion

While ugly, this approach can allow pgRouting to respect point objects with costs.

There may already be a better way to implement point costs, but this is how I did it.

It may or may not help someone out there.

Download Tables, Data, Scripts for This Article