[postgis-users] Minimum Spanning tree from MultiLineString
Bruce Rindahl
bruce.rindahl at gmail.com
Thu Jul 1 09:42:49 PDT 2021
Perfect!
Thanks
On Wed, Jun 30, 2021 at 11:32 PM Florian Nadler <florian.nadler at cybertec.at>
wrote:
> With reference to the parameter issue, try to utilize "using":
>
> EXECUTE 'create table mst_tree as (WITH t1 as (SELECT (ST_Dump($1)).geom AS geom ), t2 as(select st_snaptogrid(geom,$2) geom from t1)select geom, st_startpoint(geom) as src, st_endpoint(geom) as dest, st_length(geom) as weight from t2)' using tree,tolerance;
>
>
> CREATE OR REPLACE FUNCTION public.mst(
> tree geometry,
> tolerance float)
> RETURNS TABLE (
> mst_geom geometry,
> mst_weight double precision )
> LANGUAGE 'plpgsql' COST 100 VOLATILE PARALLEL UNSAFE ROWS 1000AS$BODY$DECLARE edge RECORD;
> e integer;
> i integer := 0 ;BEGIN EXECUTE 'drop table if exists mst_edges';
> EXECUTE 'drop table if exists mst_visited';
>
>
> EXECUTE 'drop table if exists mst_tree';
> EXECUTE 'create table mst_tree as (WITH t1 as (SELECT (ST_Dump($1)).geom AS geom ), t2 as(select st_snaptogrid(geom,$2) geom from t1) select geom, st_startpoint(geom) as src, st_endpoint(geom) as dest, st_length(geom) as weight from t2)' using tree,tolerance;
>
> EXECUTE 'create table mst_visited as (select src as node, 1 as branch from mst_tree limit 0)';
> EXECUTE 'create table mst_edges as (select geom, weight, src, dest, 1 as mst from mst_tree limit 0)';
>
> with nodes as (
> select src as node from mst_edges
> union all select dest as node from mst_edges
> ),
> d_nodes as (
> select distinct node from nodes)
> SELECT COUNT(node)
> from d_nodes into e;
>
> for edge in execute 'SELECT geom, weight, src, dest FROM mst_tree ORDER by weight asc' LOOP EXIT WHEN e = 1;
>
> WITH next_nodes as (
> select edge.src as node union all select edge.dest as node ), -- put nodes in a single column new_nodes as (
> select node from next_nodes where node not in (select node from mst_visited)
> ), -- list of nodes not yet visited existing_nodes as (
> select next_nodes.node, branch
> from next_nodes,
> mst_visited
> where next_nodes.node = mst_visited.node
> ), -- list of nodes already visited cycle_test as (
> select case when count(node) = 2 and (max(branch) = min(branch))
> then 2 -- forms a loop, do not add to MST else 1 end as mst, -- not a loop, include edge in MST max(branch) as max_branch,
> min(branch) as min_branch from existing_nodes ),
>
> v (nodes) as (
> insert into mst_visited (node, branch)
> (select node,
> case when (select count(node) from new_nodes) = 2 then e --next_branch.n_b -- two new nodes, new branch else (select branch from existing_nodes) -- add node to existing branch end as branch from new_nodes--,next_branch )
> returning * ),
>
> v_c(node) as (
> update mst_visited set branch = d.min_branch from (select * from cycle_test) d where mst_visited.branch = d.max_branch returning * ), -- if two branches, combine branches (no effect if on same branch) This is the critical step to loop detection-- update edge table, set mst to 1 or 2 u (node) as (
> insert into mst_edges (geom, weight, src, dest, mst)
> (select edge.geom, edge.weight, edge.src, edge.dest, mst from cycle_test where mst = 1)
> returning * )
>
> select mst from cycle_test into i;
>
> IF (i = 1) then e = e - 1;
> END IF;
> END LOOP;
>
> EXECUTE 'drop table mst_visited';
> RETURN QUERY SELECT geom, weight from mst_edges;
> END;$BODY$;
> select mst(st_geometryfromtext(
> 'MULTILINESTRING((12.927 149.924,50.261 222.857),(12.927 149.924,122.094 70.04),(12.927 149.924,49.247 169.665),(49.26 169.665,178.334 165.312),(50.01 222.79,49.26 169.665),(64.208 197.562,173.254 269.374),(64.355 197.598,50.01 222.455),(64.355 197.598,178.499 165.263),(64.602 197.598,49.26 169.665),(122.093 70.04,49.26 169.665),(173.254 269.374,50.01 222.79),(173.254 269.374,232.314 286.445),(173.427 269.374,242.527 269.556),(178.427 165.263,242.527 269.556),(178.499 165.263,122.093 70.04),(178.499 165.263,173.254 269.374),(232.314 286.445,242.527 269.556),(232.314 286.437,285.906 259.373),(242.527 269.556,285.906 259.236),(285.906 259.373,122.026 70.04),(285.906 259.373,178.427 165.263))'),
> 10);
>
> Am 01.07.2021 um 00:39 schrieb Bruce Rindahl:
>
> Warning - Long post
>
> I have created a function that finds the minimum spanning tree from a
> MultiLineString per the description at:
>
> https://docs.google.com/presentation/d/1iqTLwt95rEBISBcPoaA31_clqQBAJlaMegHNfNYHuuI/edit#slide=id.g6c6fcfd43d_0_85
> I have it working in two parts but I can't merge them together into one
> function.
>
> First the data. I created a table wiki with just a geometry column. Then
> I grabbed the SVG coordinates from the Wiki link from above, tweaked it
> and got the following MultiLineString:
>
> MULTILINESTRING((12.927 149.924,50.261 222.857),(12.927 149.924,122.094
> 70.04),(12.927 149.924,49.247 169.665),(49.26 169.665,178.334
> 165.312),(50.01 222.79,49.26 169.665),(64.208 197.562,173.254
> 269.374),(64.355 197.598,50.01 222.455),(64.355 197.598,178.499
> 165.263),(64.602 197.598,49.26 169.665),(122.093 70.04,49.26
> 169.665),(173.254 269.374,50.01 222.79),(173.254 269.374,232.314
> 286.445),(173.427 269.374,242.527 269.556),(178.427 165.263,242.527
> 269.556),(178.499 165.263,122.093 70.04),(178.499 165.263,173.254
> 269.374),(232.314 286.445,242.527 269.556),(232.314 286.437,285.906
> 259.373),(242.527 269.556,285.906 259.236),(285.906 259.373,122.026
> 70.04),(285.906 259.373,178.427 165.263))
>
> I then inserted it into the table wiki. Note the vertices do not
> perfectly line up..
>
> Step one - I need to create a table that defines the geometry,starting and
> ending points of each line segment and the line length. This is easy:
>
> create table mst_tree as (
> WITH t1 as (SELECT (ST_Dump(geom)).geom AS geom from wiki),
> t2 as(select st_snaptogrid(geom,10) geom from t1)
> select geom, st_startpoint(geom) as src, st_endpoint(geom) as dest,
> st_length(geom) as weight from t2
> )
>
> Note the value of 10 in the st_snaptogrid - this assures all the close
> nodes snap to each other and will need to be a parameter to the function.
>
> Next is the function that finds the Minimum Spanning Tree. I used the
> Kruskal method - see:
>
> https://www.geeksforgeeks.org/kruskals-minimum-spanning-tree-algorithm-greedy-algo-2/
>
> The full function is:
> ---------------------------------------------------------------
>
> -- FUNCTION: public.mst()
>
> -- DROP FUNCTION public.mst();
>
> CREATE OR REPLACE FUNCTION public.mst(
> -- tree geometry,
> -- tolerance float
> )
> RETURNS TABLE(mst_geom geometry, mst_weight double precision)
> LANGUAGE 'plpgsql'
> COST 100
> VOLATILE PARALLEL UNSAFE
> ROWS 1000
>
> AS $BODY$
> DECLARE
> edge RECORD;
> e integer;
> i integer := 0 ;
> BEGIN
> EXECUTE 'drop table if exists mst_edges';
> EXECUTE 'drop table if exists mst_visited';
>
> -- I am stuck here
> --EXECUTE 'drop table if exists mst_tree';
> --EXECUTE 'create table mst_tree as (WITH t1 as (SELECT (ST_Dump(tree)) AS
> geom ),
> --t2 as(select st_snaptogrid(geom,tolerance) geom from t1)
> --select geom, st_startpoint(geom) as src, st_endpoint(geom) as dest,
> st_length(geom) as weight from t2)';
>
> EXECUTE 'create table mst_visited as (select src as node, 1 as branch
> from mst_tree limit 0)';
> EXECUTE 'create table mst_edges as (select geom, weight, src, dest, 1
> as mst from mst_tree limit 0)';
>
> with nodes as (
> select src as node from mst_edges
> union all
> select dest as node from mst_edges
> ),
> d_nodes as (
> select distinct node from nodes)
> SELECT COUNT(node) from d_nodes into e;
>
> for edge in execute 'SELECT geom, weight, src, dest FROM mst_tree ORDER by
> weight asc'
> LOOP
> EXIT WHEN e = 1;
>
> WITH next_nodes as (
> select edge.src as node
> union all
> select edge.dest as node
> ), -- put nodes in a single column
>
> new_nodes as (
> select node from next_nodes where node not in (select node from
> mst_visited)
> ), -- list of nodes not yet visited
>
> existing_nodes as (
> select next_nodes.node, branch from next_nodes, mst_visited
> where next_nodes.node = mst_visited.node
> ), -- list of nodes already visited
>
> cycle_test as (
> select
> case
> when count(node) = 2 and (max(branch) = min(branch)) then 2 -- forms a
> loop, do not add to MST
> else 1 end as mst, -- not a loop, include edge in MST
> max(branch) as max_branch, min(branch) as min_branch
> from existing_nodes
> ),
>
> v (nodes) as (
> insert into mst_visited (node,branch)
> (select node,
> case
> when (select count(node) from new_nodes) = 2 then e --next_branch.n_b --
> two new nodes, new branch
> else (select branch from existing_nodes) -- add node to existing branch
> end
> as branch
> from new_nodes--,next_branch
> )
> returning *
> ),
>
> v_c(node) as (
> update mst_visited set branch = d.min_branch
> from (select * from cycle_test) d
> where mst_visited.branch = d.max_branch
> returning *
> ), -- if two branches, combine branches (no effect if on same branch)
> This is the critical step to loop detection
>
> -- update edge table, set mst to 1 or 2
> u (node) as (
> insert into mst_edges (geom,weight,src,dest,mst)
> (select edge.geom,edge.weight, edge.src, edge.dest, mst from cycle_test
> where mst = 1)
> returning *
> )
>
> select mst from cycle_test into i;
>
> IF (i = 1) then
> e = e - 1;
> END IF;
> END LOOP;
>
> EXECUTE 'drop table mst_visited';
> RETURN QUERY SELECT geom,weight from mst_edges;
>
> END;
> $BODY$;
>
> ALTER FUNCTION public.mst()
> OWNER TO postgres;
> -----------------------------------------------
>
> If you create the table mst_tree from the first part and then run:
>
> SELECT * from mst()
>
> you will get a table of each segment in the MST and the associated weight
> (length) . If you want to get the result in a MultiLine segment:
>
> SELECT ST_Collect(ARRAY(SELECT mst_geom FROM mst6()));
>
> Now I am stuck. How do I pass the geometry and tolerance into the
> function If I uncomment the lines I get the error:
>
> ERROR: column "geom" does not exist LINE 1: ...te table mst_tree as (WITH
> t1 as (SELECT (ST_Dump(geom)).geo... QUERY: create table mst_tree as
> (WITH t1 as (SELECT (ST_Dump(geom)).geom AS geom), t2 as(select
> st_snaptogrid(geom,$2) geom from t1) select geom, st_startpoint(geom) as
> src, st_endpoint(geom) as dest, st_length(geom) as weight from t2)
> Any help will be appreciated
>
>
>
>
>
>
>
> _______________________________________________
> postgis-users mailing listpostgis-users at lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/postgis-users
>
> --
> CYBERTEC PostgreSQL International GmbH
> Gröhrmühlgasse 26, A-2700 Wiener Neustadt
> Web: https://www.cybertec-postgresql.com
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210701/23c4c585/attachment.html>
More information about the postgis-users
mailing list