4

I am trying to write an SQL that calculates all points at a given distance from a point within a network (most likely the wrong term) where the network is stored in a PostgreSQL/PostGIS database as geography objects.

To try to simplify the problem I have created a network of MultilineStrings (as I get multilinestring from the database):

"MULTILINESTRING((-3395658.95218449 9999999.63159823,-3395629.87976969 9999999.12378312))"
"MULTILINESTRING((-3395640.58204213 9999999.31072238,-3395605.60319134 10000047.4596158))"
"MULTILINESTRING((-3395629.87976969 9999999.12378312,-3395601.369858 9999960.94086569))"
"MULTILINESTRING((-3395654.02194139 9999999.54548045,-3395654.02194139 9999977.87419904))"
"MULTILINESTRING((-3395654.02194139 9999977.87419904,-3395568.2969413 9999977.87419904))"
"MULTILINESTRING((-3395591.32900369 9999977.87419904,-3395591.58027466 10000017.8262824))"
"MULTILINESTRING((-3395591.46379143 9999999.30544906,-3395591.46379143 9999999.30544906,-3395618.83235802 9999999.30544906))"
"MULTILINESTRING((-3395608.51010251 9999999.30544906,-3395608.51360801 10000016.5033657))"
"MULTILINESTRING((-3395591.46379143 9999999.30544906,-3395567.93313922 9999999.33852198))"
"MULTILINESTRING((-3395582.44960782 9999999.31830406,-3395582.22063923 9999985.58018863))"
"MULTILINESTRING((-3395658.95218449 9999999.63159823,-3395684.61438934 9999974.46768862))"
"MULTILINESTRING((-3395637.01606056 9999977.94918962,-3395637.25397262 9999992.45935531))"
"MULTILINESTRING((-3395625.67040625 10000019.9754834,-3395672.70813933 10000020.5051887))"

with a starting point @ -3395649.31,9999999.42 trying to find distances 50m away.

This should give a layout like below, where the green point is the starting point and the red points are the ones I would like to calculate: enter image description here

Note: if while following a path and I hit a dead-end I would do a 180 and keep going till I reach the distance and at intersection split and follow each new path.

i.e: enter image description here

Some errors in the data I have found so far, Not all lines intersect where they should be an intersection, so I also need to work out a way to extend lines to make them intersect when following a path (normally the end/start of the line in within 500mm).

I have converted my data set to topology and noded with pgrouting (@geozelot pointed me to pgrouting), what I'm still missing is how to add a starting point to the existing noded data to do the calculations from?

enter image description here

SQL for noded data set if it helps:

CREATE TABLE public.gis_test_noded (
    id bigint NOT NULL,
    old_id integer,
    sub_id integer,
    source bigint,
    target bigint,
    cost double precision DEFAULT 0,
    geom public.geometry(LineString,32755)
);

COPY public.gis_test_noded (id, old_id, sub_id, source, target, cost, geom) FROM stdin;
1   1   1   1   2   4.930995161664557   0102000020F37F0000020000006E2EE17925E849C17E0D36F4CF126341BBF9CE0223E849C1699374F1CF126341
2   1   2   2   3   13.441949393534069  0102000020F37F000002000000BBF9CE0223E849C1699374F1CF126341455B804A1CE849C11070F1E9CF126341
3   1   3   3   4   10.703904974912364  0102000020F37F000002000000455B804A1CE849C11070F1E9CF1263410B4B9CF016E849C10708F6E3CF126341
4   2   1   3   22  25.48298626673298   0102000020F37F000002000000465B804A1CE849C11070F1E9CF12634111075FCD14E849C123B6AE7DD2126341
5   2   2   22  17  34.03033945258891   0102000020F37F00000200000011075FCD14E849C123B6AE7DD2126341B45F35CD0AE849C1272CB5EED5126341
6   3   1   4   5   26.519541477588916  0102000020F37F0000020000000B4B9CF016E849C10708F6E3CF1263413189B7010FE849C14470F93BCD126341
7   3   2   5   6   21.132848264454513  0102000020F37F0000020000003189B7010FE849C14470F93BCD126341CD8157AF08E849C15D921B1ECB126341
8   5   1   7   20  17.005880832206458  0102000020F37F000002000000BBF9CE0223E849C14470F93BCD126341BC450E821AE849C14470F93BCD126341
9   5   2   20  5   23.002646987792104  0102000020F37F000002000000BC450E821AE849C14470F93BCD1263413189B7010FE849C14470F93BCD126341
10  5   3   5   8   22.6844098768197    0102000020F37F0000020000003189B7010FE849C14470F93BCD12634106CB1CAA03E849C14470F93BCD126341
11  5   4   8   9   23.0320623931475    0102000020F37F00000200000006CB1CAA03E849C14470F93BCD1263412D2C0226F8E749C14470F93BCD126341
12  6   1   8   10  21.431673879337882  0102000020F37F00000200000006CB1CAA03E849C14470F93BCD12634181845DBB03E849C11D3DC6E9CF126341
13  6   3   10  11  18.521199648420286  0102000020F37F00000200000081845DBB03E849C11D3DC6E9CF126341A67046CA03E849C1D2E7703AD2126341
14  7   1   10  12  17.046311081852764  0102000020F37F00000200000081845DBB03E849C11D3DC6E9CF126341040A4B410CE849C11D3DC6E9CF126341
15  7   2   12  13  10.322255508508533  0102000020F37F000002000000040A4B410CE849C11D3DC6E9CF12634127B58A6A11E849C11D3DC6E9CF126341
16  9   1   10  14  9.014195599909053   0102000020F37F00000200000081845DBB03E849C11D3DC6E9CF126341DFA58C39FFE749C153072EEACF126341
17  9   2   14  15  14.516479855393287  0102000020F37F000002000000DFA58C39FFE749C153072EEACF1263411B1B71F7F7E749C10C2CD5EACF126341
19  10  2   14  16  13.739838008188064  0102000020F37F000002000000DFA58C39FFE749C153072EEACF12634108E83D1CFFE749C1C0E79032CE126341
20  4   1   2   7   21.671281406655908  0102000020F37F000002000000BBF9CE0223E849C1699374F1CF126341BBF9CE0223E849C14470F93BCD126341
21  8   1   12  18  17.197917041939753  0102000020F37F000002000000040A4B410CE849C11D3DC6E9CF12634143E8BD410CE849C17C921B10D2126341
22  11  1   1   19  35.94121735091404   0102000020F37F0000020000006E2EE17925E849C17E0D36F4CF126341564FA44E32E849C1204EF7CECC126341
23  12  1   20  21  14.51211598932906   0102000020F37F000002000000BC450E821AE849C1E7C25F3ECD126341C92C82A01AE849C1E609B30ECF126341
24  13  1   22  23  47.04071556131955   0102000020F37F00000200000041DFCFD514E849C1F928377FD21263413B4FA45A2CE849C16D812A90D2126341
\.
4
  • 3
    Graph theory algorithms for routing applications are implemented in pgRouting, a PostgreSQL extension utilizing the PostGIS spatial types and functions.
    – geozelot
    May 4, 2021 at 6:05
  • Thanks @geozelot, I have read up a little on prgrouting and have converted the test data to the formats it needs, I believe I should be using pgr_drivingDistance but I'm struggling in starting from my startpoint and node the closest node, is the a way to add the startnode in to the noded dataset?
    – Yendor
    May 6, 2021 at 4:38
  • Look into pgr_withPointsDD of the still proposed, but supported withPoints - Family of functions - note that, in order to refer to the pids of the points_sql as start/end vertices, you'd have to prefix the actual values with - (minus) in the function call! See gis.stackexchange.com/questions/377630/…
    – geozelot
    May 6, 2021 at 10:04
  • ...and gis.stackexchange.com/questions/379120/… where you find a link to a post with yet some more helpful links to make use of an efficient (K)NN query to get the required edge_id and frac parameters. And after you made it work, self-answer you question!
    – geozelot
    May 6, 2021 at 10:11

1 Answer 1

2

With using pgr_withPointDD and from: https://stackoverflow.com/questions/51572251/how-to-get-a-road-path-which-can-be-travelled-in-10-minutes-from-a-location I came up with:

WITH
  dd AS (
    SELECT pg.node,
           pg.edge,
           pg.cost
    FROM pgr_withPointsDD(
        'SELECT id, source, target, cost FROM gis_test_noded ORDER BY id',
        'WITH 
                start_vertex AS (
                        SELECT ST_SetSRID(ST_MakePoint(-3395649.31, 9999999.42), 32755) AS pt
                ),
                nearest_line AS (
                        SELECT id, ST_Distance(geom, start_vertex.pt)
                        FROM gis_test_noded, start_vertex
                        ORDER BY 2 ASC LIMIT 1
                ),
                new_vertex AS (
                        SELECT nearest_line.id AS edge_id, ST_LineLocatePoint(gis_test_noded.geom, start_vertex.pt) as fraction
                        FROM gis_test_noded, start_vertex, nearest_line
                        WHERE gis_test_noded.id = nearest_line.id
                )
        SELECT edge_id, fraction FROM new_vertex',
        -1, 50,
        directed := false) AS pg
  ),
  dd_edgs AS (
    SELECT edg.id,
           edg.geom
    FROM gis_test_noded AS edg
    JOIN dd AS d1
      ON edg.source = d1.node
    JOIN dd AS d2
      ON edg.target = d2.node
  ),
  dd_ext AS (
    SELECT edg.id,
             CASE
               WHEN dd.node = edg.source
                  THEN ST_LineInterpolatePoint(edg.geom, (50 - dd.cost) / 50)
                  ELSE ST_LineInterpolatePoint(edg.geom, 1 - ((50 - dd.cost) / 50))
             END AS geom
    FROM dd
    JOIN gis_test_noded AS edg
      ON dd.node IN (edg.source, edg.target) AND edg.id NOT IN (SELECT id FROM dd_edgs)
  )

SELECT id,
       geom
FROM dd_ext;

which gives me:

"id" "geom"
5 "0101000020F37F0000603305E60FE849C18543AF2DD4126341"
24 "0101000020F37F0000283E115E20E849C17AEA8687D2126341"
23 "0101000020F37F0000B2A626961AE849C1012FC670CE126341"
7 "0101000020F37F00008990B6090CE849C19541823DCC126341"
10 "0101000020F37F0000F95B28AE09E849C14470F93BCD126341"

enter image description here

Where -1 green dot is the starting point, yellow dots are the calculated endpoints and red dots are where I would have expected to see a yellow. But for now this answers my original question of what type of SQL I would need with the added help of using pgrouting.

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.