5

I'd like to calculate the distance between one and multiple other points along a network of LineStrings. I'm aware of this question which calculates the distance on a substring of a LineString, but not over multiple ones.

Q How can I calculate the distances along the network from the Point SITE1 to the two Points WWTP1 and WWTP2 in the example data below? Is this something for pgRouting? Note, I'm not looking for the Euclidean distance.

Data

DROP TABLE IF EXISTS river;
CREATE TABLE river (gid varchar, geom geometry);
INSERT INTO river VALUES
    ('ID1', 'LINESTRING(1 1, 2 2, 2 3, 4 3)'),
    ('ID2', 'LINESTRING(4 3, 5 6, 6 8)'),
    ('ID3', 'LINESTRING(6 8, 8 4, 9 5)'),
    ('ID4', 'LINESTRING(2 3, 1 4)'),
    ('ID5', 'LINESTRING(1 4, 1 5, 0 4, -1 1)'),
    ('ID6', 'LINESTRING(-1 1, -3 2)');
    
DROP TABLE IF EXISTS site;
CREATE TABLE site (gid varchar, geom geometry);
INSERT INTO site VALUES
    ('SITE1', 'POINT(1 1)');

DROP TABLE IF EXISTS wwtp;
CREATE TABLE wwtp (gid varchar, geom geometry);
INSERT INTO wwtp VALUES
    ('WWTP1', 'POINT(0 4)'),
    ('WWTP2', 'POINT(8 4)');

2 Answers 2

5

As you already pointed out, this indeed can be solved with pgRouting and will have to be added to your database as an extension. You will have to do some preparations with the data to get it going.

First of all you will have to break the linestrings into individual segments in order to build a sound topology:

-- https://blog.cleverelephant.ca/2015/02/breaking-linestring-into-segments.html
CREATE TABLE river_segments AS
WITH 
dumps AS ( 
  SELECT gid, ST_DumpPoints(geom) AS pt FROM river
), 
pts AS (
  SELECT gid, (pt).geom, (pt).path[1] AS vert FROM dumps
) 
SELECT a.gid as source_id, row_number() over() as gid, ST_MakeLine(ARRAY[a.geom, b.geom]) AS geom
FROM pts a, pts b 
WHERE a.gid = b.gid AND a.vert = b.vert-1 AND b.vert > 1;

pgRouting requires source and target columns to compute the toplogy.

ALTER TABLE river_segments ADD COLUMN IF NOT EXISTS source integer;
ALTER TABLE river_segments ADD COLUMN IF NOT EXISTS target integer;

SELECT pgr_createTopology('river_segments', 0.0001, 'geom', 'gid');

The function pgr_createTopology will add a new table named river_segments_vertices_pgr holding your graph nodes. Now you have everything you need to compute shortest paths between 2 points. The following common table expression selects a start and your destination point. These don't have to sit directly on the topology. When we invoke pgr_dijkstra we tell PostGIS to find the nearest node as we limit it to 1. We do this for the source and destination point. Afterwards we join the shortest path results again with the river_segments table to obtain the geometries of the route. We stitch these together with ST_Collect and voila.

WITH start_point AS (
    SELECT geom from site
),
destination_point AS (
    SELECT geom from wwtp
    WHERE gid = 'WWTP1'
),
dijkstra AS (
    SELECT * FROM pgr_dijkstra(
        'SELECT gid as id, source, target, ST_Length(geom) AS cost FROM river_segments',
        -- finding the source node spatially
        (SELECT a.id FROM river_segments_vertices_pgr a, start_point b
            ORDER BY a.the_geom <-> b.geom LIMIT 1),
        -- finding the target node spatially
        (SELECT a.id FROM river_segments_vertices_pgr a, destination_point b
            ORDER BY a.the_geom <-> b.geom LIMIT 1),
        false
    )
),
route_geom AS (
    SELECT dijkstra.seq, dijkstra.cost,
    CASE
        WHEN dijkstra.node = river_segments.source THEN geom
        ELSE ST_Reverse(geom)
    END AS route_geom
    FROM dijkstra JOIN river_segments
    ON (edge = gid) ORDER BY seq
)
SELECT ST_AsText(ST_Collect(route_geom)) FROM route_geom;

routes with pgrouting

2
  • 1
    This works for me, But how to implement 1 to many using this same data and query. what changes are to be made to get path to all destination points at once ? I need each destination path in a different row ?
    – George V
    Jan 17, 2022 at 2:17
  • 1
    @timothy-dalton I accepted my answer, since it leads to my intended result. +1 and thanks though, your code helped me a lot!
    – andschar
    Mar 18, 2022 at 14:30
1

Here's an answer if you want to solve the One-to-Many situation, as I originally intended. It's a minor modification of the last code chunk of @timothy-dalton's (+1) answer. The most important part is that you pass the target id's as an array SELECT ARRAY_AGG(a.id).

WITH
start_point AS (
    SELECT geom
    FROM site
),
destination_point AS (
    SELECT gid, geom
    FROM wwtp
    WHERE gid IN ('WWTP1', 'WWTP2')
),
dijkstra AS (
    SELECT * FROM pgr_dijkstra(
        'SELECT gid as id, source, target, ST_Length(geom) AS cost FROM river_segments',
        -- finding the source node spatially
        (SELECT a.id
         FROM river_segments_vertices_pgr a, start_point b
         ORDER BY a.the_geom <-> b.geom LIMIT 1),
        -- finding all targets node spatially
        (SELECT ARRAY_AGG(a.id)
         FROM river_segments_vertices_pgr a, destination_point b
         WHERE ST_Intersects(a.the_geom, b.geom)),
        false
    )
),
route_geom AS (
    SELECT
        dijkstra.seq,
        dijkstra.cost,
        dijkstra.end_vid,
    CASE
        WHEN dijkstra.node = river_segments.source THEN geom
        ELSE ST_Reverse(geom)
    END AS route_geom
    FROM dijkstra
    LEFT JOIN river_segments ON (edge = gid)
    ORDER BY seq
),
dist_geom AS (
    SELECT
        rou.end_vid,
        ST_AsText(ST_Collect(route_geom)) wkt,
        ST_Length(ST_Collect(route_geom)) length,
        ST_Collect(route_geom) geom
    FROM route_geom rou
    GROUP BY rou.end_vid
)
SELECT wwt.gid, dis.*
FROM dist_geom dis
LEFT JOIN wwtp wwt ON ST_Intersects(dis.geom, wwt.geom)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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