Geographic Information Systems Asked by Yendor on July 29, 2021
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:
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.
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?
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
.
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" |
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.
Correct answer by Yendor on July 29, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP