pgr_withPointsCost returns incorrect table values

Question

I have been struggling with getting pgr_withPointsCost to output a correct cost table for months. I have created a routable network, ensured the network is connected, created a poi table with all the required fields, but it only outputs a single output whose distance isn't correct. Everything is in EPSG 3400.

For the poi table (migration.cntrlpts_move), I have moved all the points to the line network, set the pid column to have a negative value as required. It should start at id = 1 (pid=-1) and end at id = 6 (pid = -6)

id	fraction	edge_id	pid
1	0	1790	-1
2	1	1792	-2
3	0.98	3319	-3
4	0.27	2244	-4
5	0.406	3523	-5
6	0	4020	-6

The output table looks like this:

start_pid	end_pid	agg_cost	id	fraction	geom	edge_id	pid
1	2	215.0217	2	1	0101000020480D000083BD86371B15144165D5F8130F375741	1792	2

The following are all the steps:

Step 1: Create Routable Network

-- ### Step 1 - Create Topology of imported streams layer ###
-- Build topology for stream CCP distance calculations
ALTER TABLE migration.streams ADD COLUMN "source" integer; -- add source column for pg routing
ALTER TABLE migration.streams ADD COLUMN "target" integer; -- add target column for pg routing
SELECT pgr_createTopology('migration.streams', 1.0, 'geom', 'id');

-- ### Step 2 - Add pgRouting networking information
ALTER TABLE migration.streams ADD COLUMN length float8; -- add length field for pgr cost
ALTER TABLE migration.streams ADD COLUMN cost float8; -- add cost column
ALTER TABLE migration.streams ADD COLUMN reverse_cost float8; -- add reverse_cost column

UPDATE migration.streams SET length = ST_Length(geom);
UPDATE migration.streams SET cost = ST_length(geom); 
UPDATE migration.streams SET reverse_cost = -1; 

-- ### Step 4 - Move release and CCP locations to stream and calculate fraction of stream CCP is located on
-- We first need to get a candidate set of maybe-closest stream, ordered by id and distance.
DROP TABLE IF EXISTS migration.control_pts_events;
CREATE TABLE migration.control_pts_events AS
WITH ordered_nearest AS (
SELECT
  ST_GeometryN(streams.geom,1) AS streams_geom,
  streams.id AS streams_id,
  pts.geom AS pts_geom,
  pts.id AS pts_id,
  pts.ccp_no,
  pts.pipeline,
  pts.vel_lo,
  pts.vel_hi,
  ST_Distance(streams.geom, pts.geom) AS distance_to_stream
FROM migration.streams streams
  JOIN migration.ccp pts                        -- *** change to shapefile name ***
  ON ST_DWithin(streams.geom, pts.geom, 500)    -- *** search distance away from stream ***
ORDER BY pts_id, distance_to_stream ASC
)
SELECT
  DISTINCT ON (pts_id)
  pts_id,
  streams_id,
  ST_LineLocatePoint(streams_geom, pts_geom) AS fraction,
  distance_to_stream
FROM ordered_nearest;


-- Table for pgr_withPointsCost operation with the fraction of the stream a point is on & stream edge_id.
DROP TABLE IF EXISTS migration.cntrlpts_move;
CREATE TABLE migration.cntrlpts_move AS
SELECT
  events.pts_id AS id, events.fraction,
  ST_LineInterpolatePoint(ST_GeometryN(streams.geom, 1), events.fraction) AS geom,
  events.streams_id AS edge_id
FROM migration.control_pts_events events
JOIN migration.streams streams
ON (streams.id = events.streams_id);

-- Update control points table to include -pid
ALTER TABLE migration.cntrlpts_move ADD COLUMN "pid" bigint;        -- add pid column for pg routing
UPDATE migration.cntrlpts_move SET pid = id;


-- Cost path from ccp to ccp
DROP TABLE IF EXISTS migration.ccp_cost;
CREATE TABLE migration.ccp_cost AS
WITH pgr_ccp AS (
SELECT r.*, p.*
FROM   pgr_withPointsCost(
         'SELECT id, 
                 source, 
                 target, 
                 length AS cost, 
                 reverse_cost
                 FROM migration.streams',
         'SELECT pid*-1 AS pid, 
                 edge_id, 
                 fraction 
                 FROM migration.cntrlpts_move WHERE pid < 0',
            (SELECT pid FROM migration.cntrlpts_move WHERE id = 1),             -- start location id number
            ARRAY(SELECT pid FROM migration.cntrlpts_move WHERE id > 1),        -- not equal to start location
            directed => true) AS r INNER JOIN migration.cntrlpts_move AS p ON
            r.end_pid = p.pid;
    ), ccp_distances AS (
    SELECT * FROM pgr_ccp 
    )
    SELECT * FROM ccp_distances ORDER BY agg_cost;

Answer 1

Update, having access to exemplary data:

While my answer is still valid for your presented query, there is another issue with your setup that is not very well documented: your Point SQL produces vertices that are considered to not actually be outside of the graph, i.e. where ST_LineLocatePoint is either 0.0 or 1.0 - and pgRoutings graph engine is rather picky about this when called with the withPoints family of functions.

The general idea here is that to switch using the pid of the closest regular vertex instead of an outside pid whenever the fraction is either 0 or 1. The sophisticated solution would include a conditional check on a prepared POI table to identify outside points that project onto regular vertices, and

• avoid them in the Points SQL with a filter
• use the respective regular pids in the from_pid/to_pid conditionally

A sufficiently precise workaround is to simply de- or increase the fraction in those cases, e.g. using my query below:

...
CASE WHEN stm._frac = 0.0 THEN 0.0000001 WHEN stm._frac = 1.0 THEN 0.9999999 ELSE stm._frac END AS fraction,
...

---

When you run

UPDATE migration.cntrlpts_move SET pid = id;

the pids are not negative - however, you assume they are, both

• in the Point SQL, where you are negating them again to be positive
  • here they have to be positive
• when passing in start_/stop_vertices, where you are not negating them
  • here they have to be negative

Seeing that the result lists only positive pids it is obvious that pgr_withPointsCost routes between regular topology vertices only - if it were routing between Point SQL vertices (only), both start_pid & end_pid in the results would be negative.

---

Notes:

• no need to actually "move" the migration.ccp geometries; the Point SQL requires only the fraction of line length at which each ccp.geom projects onto its nearest streams.geom
• use a more streamlined (and highly efficient) (K)NN query to find the closest streams.geom to each ccp.geom - some references:
  • within an UPDATE:
    How to UPDATE with LATERAL Nearest-Neighbour query?
  • generally (and also comparing to your DISTINCT ON approach):
    How to snap points to line
• be extremely cautious about using MULTI geometries - since you are using ST_GeometryN I assume this is the case - as their network topology rules are undefined

---

Bonus:

You could simplify your setup significantly after Step 1:

-- ### Step 2 - Create cost matrix
SELECT
  cmx.*,
  pts.*
FROM
  pgr_withPointsCost(
    -- Edge SQL - populating 'cost'/'reverse_cost' on the fly
    $$
    SELECT
      id,
      source,
      target,
      ST_Length(geom) AS cost,
      -1 AS reverse_cost
    FROM
      migration.streams
    $$,
    -- Point SQL - populating 'fraction' on the fly (via (K)NN sub-query)
    $$
    SELECT
      pts.id AS pid,
      stm.id AS edge_id,
      stm._frac AS fraction
    FROM
      migration.ccp AS pts
      CROSS JOIN LATERAL (
        SELECT
          _stm.id,
          ST_LineLocatePoint(_stm.geom, pts.geom) AS _frac
        FROM
          migration.streams AS _stm
        ORDER BY
          _stm.geom <-> pts.geom
        LIMIT
          1
      ) AS stm
    $$,
    -- 'start_pid' - negating on the fly
    (
      SELECT
        -id
      FROM
        migration.ccp
      WHERE
        id = 1
    ),
    -- 'end_pids' - negating on the fly
    (
      SELECT
        ARRAY_AGG(-id)
      FROM
        migration.ccp
      WHERE
        id > 1
    )
  ) AS cmx
  JOIN
  migration.ccp AS pts
    ON -(cmx.end_pid) = pts.id  -- -(-n) = n
;