Top 5 Recent Articles
ARTICLES CATEGORIES
- Algorithms (22)
- All (399)
- Biography (1)
- Blog (44)
- Business Requirements (1)
- Commentary (1)
- Conversion (2)
- Customers (2)
- Data Models (1)
- Education (2)
- GeoRaptor (13)
- GPS (1)
- Image Processing (2)
- Import Export (8)
- Licensing (2)
- LiDAR (1)
- Linear Referencing (4)
- Manifold GIS (3)
- Mapping (1)
- MySQL Spatial (7)
- Networking and Routing (including Optimization) (5)
- Open Source (18)
- Oracle Spatial and Locator (194)
- Partitioning (1)
- PostGIS (36)
- Projections (1)
- Published Articles (1)
- qGIS (1)
- Recommendations (1)
- Services (1)
- Software Change Log (1)
- Source Code (37)
- Space Curves (9)
- Spatial Database Functions (109)
- Spatial DB comparison (1)
- Spatial XML Processing (11)
- SQL Server Spatial (92)
- Standards (3)
- Stored Procedure (17)
- Tessellation or Gridding (10)
- Tools (2)
- Topological Relationships (1)
- Training (2)
- Triangulation (2)
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.
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:
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.
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.