Create Topology With Z Levels at Intersections
I have gotten some requests for help creating topologies that have support for z-levels at intersection. The best way to understand this is to think of an overpass and underpass. At the point of intersection there is a common node but we obviously don't want to route vehicles off the overpass 30 feet down to the underpass. Typically this is managed by assigning a z-level to the street segment. Consider this crude graphic:
E(0)
/
----B(1)-----
/ / \
A(0)---/ B(0) \-----C(0)
/
/
D(0)
A-B-C is a highway overpass B-B-E is a highway underpass
(n) is the zlevel
So B(0) and B(1) have the same x-y position but they are not connected together.
You can make a fairly simple clone to pgr_createtopology() to take the zlevel*factor as the z value of a point with the x-y values and check for 3D distance between the nodes. This will result in the B(0) and B(1) getting assigned different node ids and the related segments will then not be connected.
The following is some old code I have used in the past for handle this. It is not guaranteed that this will work for anyone, but I hope it will give you a starting place to get something working. Hopefully we will find some time to integrate this into v2.1.0 whenever that happens.
drop index source_idx;
drop index target_idx;
update st set source=null, target=null;
delete from geometry_columns where f_table_name='vertices_tmp';
vacuum analyze;
CREATE OR REPLACE FUNCTION point_to_id(point geometry,
tolerance double precision)
RETURNS INT AS
$$
DECLARE
row record;
point_id int;
BEGIN
LOOP
-- TODO: use && and index
SELECT INTO row id, the_geom FROM vertices_tmp WHERE
expand(point, tolerance) && the_geom and
distance(the_geom, point) < tolerance;
point_id := row.id;
IF NOT FOUND THEN
INSERT INTO vertices_tmp (the_geom) VALUES (point);
ELSE
EXIT;
END IF;
END LOOP;
RETURN point_id;
END;
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
CREATE OR REPLACE FUNCTION point_to_id_z(point geometry, zlev float,
tolerance double precision)
RETURNS INT AS
$$
DECLARE
row record;
point_id int;
pnt geometry;
BEGIN
pnt := st_translate(
st_force_3d(point), 0.0, 0.0, coalesce(zlev::float, 0.0));
LOOP
SELECT INTO row id, the_geom FROM vertices_tmp WHERE
expand(pnt, tolerance) && the_geom and
st_length3d(makeline(the_geom, pnt)) < tolerance;
point_id := row.id;
IF NOT FOUND THEN
INSERT INTO vertices_tmp (the_geom) VALUES (pnt);
ELSE
EXIT;
END IF;
END LOOP;
RETURN point_id;
END;
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
DROP FUNCTION add_vertices_geometry(geom_table varchar);
/*
CREATE OR REPLACE FUNCTION add_vertices_geometry(geom_table varchar)
RETURNS VOID AS
$$
DECLARE
vertices_table varchar := quote_ident(geom_table) || '_vertices';
BEGIN
BEGIN
EXECUTE 'SELECT addGeometryColumn(''' ||
quote_ident(vertices_table) ||
''', ''the_geom'', -1, ''POINT'', 2)';
EXCEPTION
WHEN DUPLICATE_COLUMN THEN
END;
EXECUTE 'UPDATE ' || quote_ident(vertices_table) ||
' SET the_geom = NULL';
EXECUTE 'UPDATE ' || quote_ident(vertices_table) ||
' SET the_geom = startPoint(geometryn(m.the_geom, 1)) FROM ' ||
quote_ident(geom_table) ||
' m where geom_id = m.source';
EXECUTE 'UPDATE ' || quote_ident(vertices_table) ||
' set the_geom = endPoint(geometryn(m.the_geom, 1)) FROM ' ||
quote_ident(geom_table) ||
' m where geom_id = m.target_id AND ' ||
quote_ident(vertices_table) ||
'.the_geom IS NULL';
RETURN;
END;
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
*/
CREATE OR REPLACE FUNCTION assign_vertex_id_zlev(geom_table varchar,
f_zlev text,
t_zlev text,
tolerance double precision,
geo_cname varchar,
gid_cname varchar)
RETURNS VARCHAR AS
$$
DECLARE
points record;
source_id int;
target_id int;
srid integer;
countids integer;
sql text;
BEGIN
BEGIN
DROP TABLE vertices_tmp;
EXCEPTION
WHEN UNDEFINED_TABLE THEN
END;
EXECUTE 'CREATE TABLE vertices_tmp (id serial)';
EXECUTE 'SELECT coalesce(srid, 4326) as srid FROM geometry_columns WHERE f_table_name='|| quote_literal(geom_table) INTO srid;
EXECUTE 'SELECT count(*) as countids FROM '|| quote_ident(geom_table) INTO countids;
EXECUTE 'SELECT addGeometryColumn(''vertices_tmp'', ''the_geom'', '||srid||', ''POINT'', 3)';
CREATE INDEX vertices_tmp_idx ON vertices_tmp USING GIST (the_geom);
FOR points IN EXECUTE 'SELECT ' || quote_ident(gid_cname)
|| ' AS id, coalesce('
|| quote_ident(f_zlev) || '::float,0.0) AS f_zlev, coalesce('
|| quote_ident(t_zlev) || '::float,0.0) AS t_zlev, '
|| ' PointN('|| quote_ident(geo_cname) ||', 1) AS source,'
|| ' PointN('|| quote_ident(geo_cname) ||', NumPoints('
|| quote_ident(geo_cname) ||')) AS target'
|| ' FROM ' || quote_ident(geom_table) || ' ORDER BY '
|| quote_ident(gid_cname)
loop
IF points.id%1000=0 THEN
RAISE NOTICE '% out of % edges processed', points.id, countids;
END IF;
source_id := point_to_id_z(setsrid(points.source, srid),
points.f_zlev, tolerance);
target_id := point_to_id_z(setsrid(points.target, srid),
points.t_zlev, tolerance);
if source_id is null or target_id is null then
raise notice 'source/target is null id=%', points.id;
raise notice 'source=%, f_zlev=%', astext(points.source), points.f_zlev;
raise notice 'target=%, t_zlev=%', astext(points.target), points.t_zlev;
end if;
sql := 'update ' || quote_ident(geom_table)
|| ' SET source = ' || source_id
|| ', target = ' || target_id
|| ' WHERE ' || quote_ident(gid_cname) || ' = '
|| points.id;
--RAISE NOTICE 'source_id=%, target_id=%, sql=%', source_id, target_id, sql;
EXECUTE sql;
END LOOP;
RETURN 'OK';
END;
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
CREATE OR REPLACE FUNCTION assign_vertex_id(geom_table varchar,
tolerance double precision,
geo_cname varchar,
gid_cname varchar)
RETURNS VARCHAR AS
$$
DECLARE
points record;
source_id int;
target_id int;
srid integer;
countids integer;
BEGIN
BEGIN
DROP TABLE vertices_tmp;
EXCEPTION
WHEN UNDEFINED_TABLE THEN
END;
EXECUTE 'CREATE TABLE vertices_tmp (id serial)';
EXECUTE 'SELECT srid FROM geometry_columns WHERE f_table_name='|| quote_ident(geom_table) INTO srid;
EXECUTE 'SELECT count(*) as countids FROM '|| quote_ident(geom_table) INTO countids;
EXECUTE 'SELECT addGeometryColumn(''vertices_tmp'', ''the_geom'', '||srid||', ''POINT'', 3)';
CREATE INDEX vertices_tmp_idx ON vertices_tmp USING GIST (the_geom);
FOR points IN EXECUTE 'SELECT ' || quote_ident(gid_cname) || ' AS id,'
|| ' PointN('|| quote_ident(geo_cname) ||', 1) AS source,'
|| ' PointN('|| quote_ident(geo_cname) ||', NumPoints('
|| quote_ident(geo_cname) ||')) AS target'
|| ' FROM ' || quote_ident(geom_table) || ' ORDER BY '
|| quote_ident(gid_cname)
loop
IF points.id%1000=0 THEN
RAISE NOTICE '% out of % edges processed', points.id, countids;
END IF;
source_id := point_to_id(setsrid(points.source, srid), tolerance);
target_id := point_to_id(setsrid(points.target, srid), tolerance);
EXECUTE 'update ' || quote_ident(geom_table)
|| ' SET source = ' || source_id
|| ', target = ' || target_id
|| ' WHERE ' || quote_ident(gid_cname) || ' = '
|| points.id;
END LOOP;
RETURN 'OK';
END;
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
begin;
--select assign_vertex_id('st', 0.000001, 'the_geom', 'gid');
select assign_vertex_id_zlev('st', 'f_zlev', 't_zlev', 0.000001, 'the_geom', 'gid');
create unique index vertices_tmp_id_idx on vertices_tmp using btree(id);
create index source_idx on st using btree(source);
create index target_idx on st using btree(target);
commit;
--vacuum analyze;
--select count(*) from (select a.id, count(*) as cnt from vertices_tmp a, st b where a.id=b.source or a.id=b.target group by a.id having count(*)<2) as foo;
alter table vertices_tmp add column cnt integer;
update vertices_tmp as a set cnt=(select count(*) from st b where a.id=b.source or a.id=b.target);
vacuum analyze vertices_tmp;
If you want to help or fund this please contact us. -Stephen Woodbridge
