Friday, 20 January 2017

postgis - ST_LineMerge to simplify road network


I am trying to simplify my road network by finding all segments that touch 2 other lines with a common attribute, join_id and street, then doing a st_linemerge. I do NOT want to merge lines where there is a legitimate intersection. The end goal is to remove any unnecessary nodes in my road network.


The red circles identify the line segments that should be merged into 1, resulting in 8 linestrings:


enter image description here


Sample data:


CREATE TABLE edges_test (
id integer,

join_id character varying(15),
street character varying(32),
the_geom geometry(MultiLineString,4326)
);

INSERT INTO edges_test
(
id, join_id, street, the_geom
)
VALUES

(152919,'5364001050000','HYLAN BOULEVARD',ST_GeomFromText('MULTILINESTRING((-74.076549037205 40.6062177712225,-74.0765673473598 40.6061216642362))', 4326)),
(163162,'5364001050000','HYLAN BOULEVARD',ST_GeomFromText('MULTILINESTRING((-74.0765673473598 40.6061216642362,-74.0765887729004 40.6060012862247))', 4326)),
(164278,'5154601070000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.0747578790902 40.6059918638443,-74.0748779616281 40.6060160205844,-74.0753290236186 40.6060929633246,-74.0754908197269 40.6061131335313))', 4326)),
(164279,'5154601070000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.0754908197269 40.6061131335313,-74.0758750018611 40.6061610259355,-74.0760098160418 40.606172376073,-74.0760455106313 40.6061753816192))', 4326)),
(164280,'5154601070000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.0760455106313 40.6061753816192,-74.076549037205 40.6062177712225))', 4326)),
(164284,'5154601060000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.0748137615045 40.6057770701806,-74.0752209528085 40.6058469663565,-74.0755409340908 40.6058951565944))', 4326)),
(164285,'5154601060000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.0755409340908 40.6058951565944,-74.0757990415962 40.6059340275554,-74.0758941551205 40.6059421283636,-74.0760898820693 40.6059587982079))', 4326)),
(164286,'5154601060000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.0760898820693 40.6059587982079,-74.0765887729004 40.6060012862247))', 4326)),
(166736,'5364001050000','HYLAN BOULEVARD',ST_GeomFromText('MULTILINESTRING((-74.0765017287402 40.606466130529,-74.0765330374503 40.6063017649758))', 4326)),
(166737,'5364001050000','HYLAN BOULEVARD',ST_GeomFromText('MULTILINESTRING((-74.0765330374503 40.6063017649758,-74.076549037205 40.6062177712225))', 4326)),

(166738,'5209101000000','STATEN ISLAND EXPWY ET 13 B WB',ST_GeomFromText('MULTILINESTRING((-74.076549037205 40.6062177712225,-74.0770963146923 40.6063461215705))', 4326)),
(166740,'5154601070000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.076549037205 40.6062177712225,-74.077116909924 40.6062540917328))', 4326)),
(166742,'5154601060000','STATEN ISLAND EXPRESSWAY',ST_GeomFromText('MULTILINESTRING((-74.0765887729004 40.6060012862247,-74.0771604232288 40.6060339913006,-74.0771785369999 40.6060348585715))', 4326)),
(166851,'5364001050000','HYLAN BOULEVARD',ST_GeomFromText('MULTILINESTRING((-74.0765887729004 40.6060012862247,-74.0766012139353 40.6059313624662))', 4326))
;

My (very) wrong query:


SELECT 
e1.join_id, -- identifier for same street segment
e1.street,

ST_LineMerge(St_union(e1.the_geom)) as the_geom,
COUNT(*)
INTO edges_output
FROM edges_test e1
JOIN edges_test e2 ON ST_Touches(e1.the_geom, e2.the_geom)
GROUP BY
e1.join_id,
e1.street
HAVING count(*) = 2;

Answer




would something like this do the job?


UPDATE1: used st_crosses instead of st_intersects and not st_equals and created an index for the geometry


UPDATE2: instead of calculating a huge multipoint, lets calculate one multipoint per join_id and let's split only the line with this join_id. Also, if a particular join_id has only one segment bring it anyways


UPDATE3: update2 had a bug -- this fixes it. morespecifically it didn't take into account that some lines would touch instead of cross.


 WITH edges_test_merged as (SELECT join_id, 
street,
ST_LineMerge(ST_Union(f.the_geom)) as singlegeom
FROM edges_test As f
GROUP BY 1,2
),


edges_intersections as( SELECT a.join_id, a.street, st_union(st_intersection(a.singlegeom, b.singlegeom)) as geom
FROM edges_test_merged a
JOIN edges_test_merged b
ON ST_crosses(a.singlegeom, b.singlegeom) or st_touches(a.singlegeom, b.singlegeom)
group by 1,2)

SELECT l.join_id as ljoinid, l.street as lstreet
,(st_dump(st_split(l.singlegeom, p.geom))).geom
FROM edges_test_merged as l

JOIN edges_intersections as p
ON l.join_id = p.join_id

In addition you can index your geometry column (don't know if you already have) by


CREATE INDEX edges_test_geom_idx ON edges_test USING GIST (the_geom);

Query cost before:


Result  (cost=21024.52..21348.80 rows=1000 width=162) (actual time=1.481..1.673 rows=8 loops=1)

Query cost after update1:



Result  (cost=55.32..327.05 rows=1000 width=162) (actual time=1.297..1.436 rows=8 loops=1)

Query cost after update2:


Result  (cost=55.36..3804.30 rows=14000 width=162) (actual time=1.773..1.821 rows=8 loops=1)

Query cost after update3:


Result  (cost=104.35..372.45 rows=1000 width=162) (actual time=2.121..2.160 rows=8 loops=1)

probably there is a better way to do this. I am curious at to how this will perform on a large dataset.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...