9

I have a data of points along the streets, I would like to turn those dots into simple coloured lines. Any pointers what this problem may be called or any algorithms that can help me solve this? Points along the street I would like to turn into lines.

I was hoping to use PostGIS functions to do this but I'm open to suggestions, this is a data from a .shp file.

Edit1: Updated the picture to demonstrate ideal solution of this problem.

Drawing the line would be purely based on the distance between those points, there's nothing else that I can use to group them by. Ideally this would be points at max specified distance along the projected line? And by projected line I mean find 1st point then next one nearest to it then project a line and check if there are any points on this line at max distance to any of those already on the line.

  • 1
    Which software are you planning to use? – ArMoraer Apr 4 '16 at 20:06
  • are you trying to turn these into sidewalks? – DPSSpatial Apr 4 '16 at 20:12
  • I was hoping to use PostGIS functions to do this but I'm open to suggestions, this is a data from a .shp file. – Mahakala Apr 4 '16 at 20:37
  • 1
    Could you show exactly which points you want to connect on your drawing or on another drawing? Is it only two points at a time? Or three? Is the distance between the points that should be connected always the same or or is it "just" below a certain threshold? – Peter Horsbøll Møller Apr 5 '16 at 6:50
  • 1
    Big thanks to both @dbaston and MarHoff, I won't have time to test your ideas until late April, i wish I could split bounty between but I will need to award this to 1 of you and dbaston gave me some queries too look into so I will accept his response. Thx everyone who took time to answer! Great community to be part of :-) – Mahakala Apr 10 '16 at 8:20
8

You can use a recursive query to explore nearest neighbor of each point starting from each detected end of lines you want to build.

Prerequisites : prepare a postgis layer with your points and another with a single Multi-linestring object containing your roads. The two layers must be on the same CRS. Here is the code for the test data-set I created, please modify it as needed. (Tested on postgres 9.2 and postgis 2.1)

WITH RECURSIVE
points as (SELECT id, st_transform((st_dump(wkb_geometry)).geom,2154) as geom, my_comment as com FROM mypoints),
roads as (SELECT st_transform(ST_union(wkb_geometry),2154) as geom from highway),

enter image description here

Here are the steps:

  1. Generate for each point the list of every neighbors and theirs distance that meet theses three criteria.

    • Distance must not exceed a user defined threshold (this will avoid linking to isolated point) enter image description here
      graph_full as (
      SELECT a.id, b.id as link_id, a.com, st_makeline(a.geom,b.geom) as geom, st_distance(a.geom,b.geom) as distance
      FROM points a
      LEFT JOIN points b ON a.id<>b.id
      WHERE st_distance(a.geom,b.geom) <= 15
      ),
      
    • Direct path must not cross a road enter image description here
      graph as (
      SELECt graph_full.*
      FROM graph_full RIGHT JOIN
      roads ON st_intersects(graph_full.geom,roads.geom) = false
      ),
      
    • Distance must not exceed a user defined ratio of the distance from the nearest neighbor (this should accommodate better to irregular digitalization than fixed distance) This part was actually too hard to implement, sticked to fixed search radius

    Let's call this table "the graph"

  2. Select end of line point by joining to the graph and keeping only point that have exactly one entry in the graph. enter image description here

    eol as (
    SELECT points.* FROM
    points  JOIN
    (SELECT id, count(*) FROM graph 
    GROUP BY id
    HAVING count(*)= 1) sel
    ON points.id = sel.id),
    

    Let's call this table "eol" (end of line)
    easy? that the reward for doing a great graph but hold-on things will go crazy at next step

  3. Set up a recursive query that will cycle from neighbors to neighbors starting from each eol enter image description here

    • Initialize the recursive query using eol table and adding a counter for the depth, an aggregator for the path, and a geometry constructor to build the lines
    • Move to next iteration by switching to nearest neighbor using the graph and checking that you never go backward using the path
    • After the iteration is finished keep only the longest path for each starting point (if your dataset include potential intersection between expect lines that part would need more conditions)
    recurse_eol (id, link_id, depth, path, start_id, geom) AS (--initialisation
    SELECT id, link_id, depth, path, start_id, geom FROM (
        SELECT eol.id, graph.link_id,1 as depth,
        ARRAY[eol.id, graph.link_id] as path,
        eol.id as start_id,
        graph.geom as geom,
        (row_number() OVER (PARTITION BY eol.id ORDER BY distance asc))=1 as test
        FROM eol JOIn graph ON eol.id = graph.id 
        ) foo
    WHERE test = true
    
    UNION ALL ---here start the recursive part
    
    SELECT id, link_id, depth, path, start_id, geom  FROM (
        SELECT graph.id, graph.link_id, r.depth+1 as depth,
        path || graph.link_id as path,
        r.start_id,
        ST_union(r.geom,graph.geom) as geom,
        (row_number() OVER (PARTITION BY r.id ORDER BY distance asc))=1 as test
        FROM recurse_eol r JOIN graph ON r.link_id = graph.id AND NOT graph.link_id = ANY(path)) foo
    WHERE test = true AND depth < 1000), --this last line is a safe guard to stop recurring after 1000 run adapt it as needed
    

    Let's call this table "recurse_eol"

  4. Keep only longest line for each start point and remove every exact duplicate path Example : paths 1,2,3,5 AND 5,3,2,1 are the same line discovered by it's two differents "end of line"

    result as (SELECT start_id, path, depth, geom FROM
    (SELECT *,
    row_number() OVER (PARTITION BY array(SELECT * FROM unnest(path) ORDER BY 1))=1 as test_duplicate,
    (max(depth) OVER (PARTITION BY start_id))=depth as test_depth
    FROM recurse_eol) foo
    WHERE  test_depth = true AND test_duplicate = true)
    
    SELECT * FROM result
  5. Manually checks remaining errors (isolated points, overlapping lines, weirdly shaped street)


Updated as promised, I still can't figure out why sometimes recursive query don't give exact same result when starting from opposite eol of a same line so some duplicate may remain in result layer as of now.

Feel free to ask I totally get that this code need more comments. Here is the full query:

WITH RECURSIVE
points as (SELECT id, st_transform((st_dump(wkb_geometry)).geom,2154) as geom, my_comment as com FROM mypoints),
roads as (SELECT st_transform(ST_union(wkb_geometry),2154) as geom from highway),

graph_full as (
    SELECT a.id, b.id as link_id, a.com, st_makeline(a.geom,b.geom) as geom, st_distance(a.geom,b.geom) as distance
    FROM points a
    LEFT JOIN points b ON a.id<>b.id
    WHERE st_distance(a.geom,b.geom) <= 15
    ),

graph as (
    SELECt graph_full.*
    FROM graph_full RIGHT JOIN
    roads ON st_intersects(graph_full.geom,roads.geom) = false
    ),

eol as (
    SELECT points.* FROM
    points  JOIN
        (SELECT id, count(*) FROM graph 
        GROUP BY id
        HAVING count(*)= 1) sel
    ON points.id = sel.id),


recurse_eol (id, link_id, depth, path, start_id, geom) AS (
    SELECT id, link_id, depth, path, start_id, geom FROM (
        SELECT eol.id, graph.link_id,1 as depth,
        ARRAY[eol.id, graph.link_id] as path,
        eol.id as start_id,
        graph.geom as geom,
        (row_number() OVER (PARTITION BY eol.id ORDER BY distance asc))=1 as test
        FROM eol JOIn graph ON eol.id = graph.id 
        ) foo
    WHERE test = true

UNION ALL
    SELECT id, link_id, depth, path, start_id, geom  FROM (
        SELECT graph.id, graph.link_id, r.depth+1 as depth,
        path || graph.link_id as path,
        r.start_id,
        ST_union(r.geom,graph.geom) as geom,
        (row_number() OVER (PARTITION BY r.id ORDER BY distance asc))=1 as test
        FROM recurse_eol r JOIN graph ON r.link_id = graph.id AND NOT graph.link_id = ANY(path)) foo
    WHERE test = true AND depth < 1000),

result as (SELECT start_id, path, depth, geom FROM
    (SELECT *,
    row_number() OVER (PARTITION BY array(SELECT * FROM unnest(path) ORDER BY 1))=1 as test_duplicate,
    (max(depth) OVER (PARTITION BY start_id))=depth as test_depth
    FROM recurse_eol) foo
WHERE  test_depth = true AND test_duplicate = true)

SELECT * FROM result

  • Hi @MarHoff, thanks for your answer I have something to go after. I wasn't expecting a full solution, just a pointers where to look for answers. I want to understand this more and I will keep digging and probably have more questions later. I need to understand your algorithm and this will take me some time anyway :) – Mahakala Apr 10 '16 at 8:30
  • Got a working script, preview here qgiscloud.com/MarHoff/test_qgiscloud_bis a small caveat for de-duplication remain... No more bounty, no more pression i guess, so i will release release when I can. This puzzle was fun though – MarHoff Apr 11 '16 at 14:02
  • thank you @MarHoff, if I could I would have split this bounty I can't see how I can award you with any points, but big thanks for looking into this and your proof. Looks genuine :) – Mahakala Apr 11 '16 at 21:15
  • Done. Thanks for the puzzle, and sorry for ranting. If other answer did it for you then it's totally ok sometime simple is best... My answer was maybe a little bit overthinking it. Though nice example of using CTE + recursive query + Windows function + postgis on a single query ;) – MarHoff Apr 12 '16 at 9:16
8
+100

As @FelixIP points out, the first step is to find the points that will make up each line. You can do this by calling ST_ClusterWithin with your maximum separation distance:

SELECT
  row_number() OVER () AS cid, 
  (ST_Dump(geom)).geom 
FROM (
  SELECT unnest(st_clusterwithin(geom, 0.05)) AS geom 
  FROM inputs) sq

Then, you'll need to use some heuristic to build a line through all of the points in each cluster. For example, if you can assume the desired lines to be Y-monotone, you can sort the points in each cluster and feed them into ST_MakeLine. Combining that all together would look like this:

SELECT 
  ST_MakeLine(geom ORDER BY ST_Y(geom)) AS geom
FROM (
  SELECT row_number() OVER () AS cid, 
  (ST_Dump(geom)).geom FROM (
    SELECT unnest(st_clusterwithin(geom, 0.05)) AS geom 
    FROM inputs) sq) ssq 
GROUP BY cid
  • Way to go but the Y-monotone (or even switch between X/Y-monotone) approach wont work well if the data-set contain curved road. Is it the case? Ordering algorithm is the hardest part of this question IMHO. – MarHoff Apr 7 '16 at 15:26
  • @MarHoff: yes curved roads will be an issue, but I'm trying to get most of the data transformed automatically and rest will need to be done manually. Or I'll keep digging into topic more to figure out solution but it may take longer than getting someone to fix remaining data. I'll need to evaluate results to be able to decide. Thx for pointing this out! – Mahakala Apr 7 '16 at 16:55
  • Statut tuned I just thought of a trick I need to check out... – MarHoff Apr 7 '16 at 17:06
  • Is there any robust way to do this that doesn't involve trying all possible orderings of points, and finding which one gives the shortest total length? – dbaston Apr 7 '16 at 17:07
  • If these sets of points always follow the roads, you project the position of the point onto the road segment (ST_Line_Locate_Point), then order the points by the result. – travis Apr 7 '16 at 18:37

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.