If you followed the previous two posts (here and here), you should now have your road network and events all digitised and synced to each other. This post will provide the main Postgres function (split_roads()) that loops through each event and splits the network, finally generating a new table, containing the network split at the points where any of the event values change.
Before that, you will need to create three helper functions. The first one (truncate_if_exists()) will truncate the temporary tables required while running the main function so we can run the function multiple times in case something goes wrong:
--DROP FUNCTION truncate_if_exists
CREATE OR REPLACE FUNCTION events.truncate_if_exists(_table text, _schema text DEFAULT NULL)
RETURNS void
LANGUAGE plpgsql AS
$func$
------------------------------------------------------------------------------------------------------------------------------------------
--
-- Truncates a table if exists and is not empty
-- From: https://stackoverflow.com/questions/63004656/truncate-if-exists-in-psql-function-and-call-function/63004824#63004824
--
------------------------------------------------------------------------------------------------------------------------------------------
DECLARE
_qual_tbl text := concat_ws('.', quote_ident(_schema), quote_ident(_table));
_row_found bool;
BEGIN
IF to_regclass(_qual_tbl) IS NOT NULL THEN -- table exists
EXECUTE 'SELECT EXISTS (SELECT FROM ' || _qual_tbl || ')'
INTO _row_found;
IF _row_found THEN -- table is not EMPTY
EXECUTE 'TRUNCATE ' || _qual_tbl;
RAISE NOTICE 'Table truncated: %. Reset sequence', _qual_tbl;
EXECUTE 'ALTER SEQUENCE ' || _qual_tbl || '_id_seq RESTART WITH 1';
ELSE -- optional!
RAISE NOTICE 'Table exists but is empty: %', _qual_tbl;
END IF;
ELSE -- optional!
--RAISE EXCEPTION 'here %',_qual_tbl ;
RAISE NOTICE 'Table not found: %', _qual_tbl;
END IF;
END
$func$;
Then, create the create_temp_tab() function to create the temporary network tables:
CREATE OR REPLACE FUNCTION events.create_temp_tab(pi_tmp_roads_tab varchar)
RETURNS VOID
LANGUAGE plpgsql
AS $function$
------------------------------------------------------------------------------------------------------
--
-- Creates the temp tables used in split. Called from events.split_roads function
--
------------------------------------------------------------------------------------------------------
DECLARE
BEGIN
EXECUTE format ('SELECT events.truncate_if_exists($1, ''public'')')
USING pi_tmp_roads_tab;
-- Create the temp route table with the same structure as
-- the od_proskathg table plus an additional field [from_route_id]
-- which will hold the original road_id value
EXECUTE format ('CREATE TABLE IF NOT EXISTS public.%I (
id bigserial NOT NULL,
slope15 Boolean,
width integer,
ortho50 Boolean,
geom public.geometry(multilinestringm, 2100) NULL,
CONSTRAINT pk_' || pi_tmp_roads_tab || ' PRIMARY KEY (id)
)', pi_tmp_roads_tab);
END;
$function$;
Then the 3rd helper function that does the insert in the temp tables:
CREATE OR REPLACE FUNCTION insert_split_segments(pi_is_first boolean, pi_roads_tab character varying, pi_event_name character varying, pi_attr_val character varying, pi_seg_id bigint, pi_route_table character varying, pi_segment geometry, pi_seg_id_src bigint)
RETURNS void
LANGUAGE plpgsql
AS $function$
------------------------------------------------------------------------------------------------------
--
-- Inserts the new road record into the temp table
--
------------------------------------------------------------------------------------------------------
DECLARE
v_attr_datatype varchar;
col_rec record;
v_sql TEXT:='';
v_sql_col TEXT:='';
v_sql_val TEXT:='';
tmp_roads_seq TEXT;
new_id bigint;
BEGIN
tmp_roads_seq:= pi_roads_tab || '_id_seq';
-- Before inserting into the temp table we need to find the data type of
-- the attribute so we do the required casting if needed.
EXECUTE 'select data_type from information_schema.columns
where table_name = ''' || pi_roads_tab ||
''' and column_name =''' || pi_event_name ||''''
INTO v_attr_datatype;
--RAISE NOTICE 'v_attr_datatype: %', v_attr_datatype;
--
IF pi_is_first THEN
-- This is the first round of splits so we need to take the attribute from the event table
IF v_attr_datatype = 'boolean' THEN
RAISE NOTICE 'inside';
EXECUTE format ('INSERT INTO %I (%I, from_route_id, from_route_id_src, geom)
VALUES (CAST($1 AS BOOLEAN), $2, $3, ST_Multi($4))',
pi_roads_tab, pi_event_name)
USING pi_attr_val, pi_seg_id, pi_seg_id_src, pi_segment;
ELSIF v_attr_datatype = 'integer' THEN
EXECUTE format ('INSERT INTO %I (%I, from_route_id, from_route_id_src, geom)
VALUES (CAST($1 AS INTEGER), $2, $3, ST_Multi($4))',
pi_roads_tab, pi_event_name)
USING pi_attr_val, pi_seg_id, pi_seg_id_src, pi_segment;
ELSE
EXECUTE format ('INSERT INTO %I (%I, from_route_id, from_route_id_src, geom)
VALUES ($1, $2, $3, ST_Multi($4))',
pi_roads_tab, pi_event_name)
USING pi_attr_val, pi_seg_id, pi_seg_id_src, pi_segment;
END IF;
ELSE -- Copy the attributes from the route table that we use for the split
v_sql_col:= 'id, ';
EXECUTE 'SELECT nextval('''|| tmp_roads_seq || ''')' INTO new_id;
v_sql_val:= new_id::varchar || ', ';
FOR col_rec IN
EXECUTE 'SELECT column_name from information_schema.columns
WHERE table_name = ''' || pi_roads_tab ||
''' and table_schema = ''public'' AND column_name NOT IN (''id'', ''from_route_id'', ''from_route_id_src'', ''geom'') order by ordinal_position'
LOOP
--RAISE NOTICE '%', col_rec.column_name;
v_sql_col:= v_sql_col || col_rec.column_name || ','΄;
--RAISE NOTICE '%', v_sql_col;
v_sql_val:= v_sql_val || chr(10) || '(select '||col_rec.column_name || ' FROM ' || pi_route_table || ' WHERE id = ' || pi_seg_id::varchar || '),';
END LOOP;
--RAISE NOTICE '%', pi_roads_tab;
--RAISE NOTICE '%', v_sql_col;
v_sql := 'INSERT INTO ' || pi_roads_tab || ' (' || v_sql_col || ' from_route_id, from_route_id_src, geom) VALUES (' || v_sql_val || pi_seg_id::varchar || ',' || pi_seg_id_src::varchar || ',' || 'ST_Multi(ST_GeomFromText(''' || ST_AsText(pi_segment) || ''', 2100)' || '))';
--RAISE NOTICE '%', v_sql;
EXECUTE v_sql;
-- Now do the update
--RAISE NOTICE 'Do the update';
IF v_attr_datatype = 'boolean' THEN
EXECUTE format ('UPDATE %I SET %I = CAST($1 AS BOOLEAN) WHERE id= $2',
pi_roads_tab, pi_event_name)
USING pi_attr_val, new_id;
ELSIF v_attr_datatype = 'integer' THEN
EXECUTE format ('UPDATE %I SET %I = CAST($1 AS INTEGER) WHERE id= $2',
pi_roads_tab, pi_event_name)
USING pi_attr_val, new_id;
ELSE
EXECUTE format ('UPDATE %I SET %I = $1 WHERE id= $2',
pi_roads_tab, pi_event_name)
USING pi_attr_val, new_id;
END IF;
END IF;
END;
$function$;
And finally, the main split_roads() functions that does the splitting:
CREATE OR REPLACE FUNCTION split_roads()
RETURNS void
LANGUAGE plpgsql
AS $function$
------------------------------------------------------------------------------------------------------
-- The function that will run on the AFTER INSERT trigger for any event
--
------------------------------------------------------------------------------------------------------
DECLARE
events text[] = ARRAY['slope15', 'road_width', 'ortho50'];
event_name text;
tmp_roads_tab TEXT;
event_rec record;
v_event_geom geometry;
v_attr_val varchar;
v_attr_datatype varchar;
v_route_table varchar = 'roads'; -- Name of the initial route table
v_output_tab varchar = 'split_roads'; -- The name of the output table
v_route_id int8;
v_route_geom geometry;
route_rec record;
v_start_offset NUMERIC;
v_end_offset NUMERIC;
v_startM NUMERIC;
v_endM NUMERIC;
cnt int4 = 1;
v_geom_element2 geometry; --The individual geometries within the GEOMETRYCOLLECTION from the first split
v_start_pnt geometry; --The start point of the event geometry to be used in the split
v_end_pnt geometry; --The end point of the event geometry to be used in the split
v_start_loc float8;
v_end_loc float8;
temp_route_rec record;
v_route_start_pnt geometry;
v_route_startLoc float8;
v_route_end_pnt geometry;
v_route_endLoc float8;
BEGIN
-- Start looping through the event tables
FOREACH event_name IN ARRAY events
LOOP
tmp_roads_tab:= 'temp'|| cnt::varchar || '_split_od_proskathg';
IF cnt = 1 THEN -- this IS the FIRST EVENT
RAISE NOTICE '------------------------- ITERATION #% (%)---------------------------', cnt, event_name;
PERFORM create_temp_tab(tmp_roads_tab);
--
-- Now loop through the events that are located on the network
-- (i.e. where road_id column is not null)
--
FOR event_rec IN
EXECUTE format ('SELECT id, %s::varchar, road_id, geom FROM %I WHERE road_id is not null', event_name, event_name)
USING event_name
LOOP
--- Assign the value of the dynamically fetched column to v_attr_val
EXECUTE format('SELECT %I FROM %I WHERE id = $1', event_name, event_name)
INTO v_attr_val
USING event_rec.id;
-- Since this is the first event, get the route geometry for each record using the original network table
EXECUTE format ('SELECT id, geom FROM %I WHERE id= $1', v_route_table)
USING event_rec.road_id
INTO v_route_id, v_route_geom;
-- Get the event's geometry
v_event_geom:= event_rec.geom;
-- We now have the road geometry and the event geometry. Start the splitting.
-- Since this is the first event we process, we know that the temp table is empty
-- and that the event geometry, lies wholly within the road geometry
--
-- Get the start and end points of the event geometry
v_start_pnt:= ST_StartPoint(ST_GeometryN(v_event_geom,1));
v_end_pnt:= ST_EndPoint(ST_GeometryN(v_event_geom,1));
-- Get start and end point locations along the route geometry
-- IMPORTANT NOTE: Always assume single part route geometries
v_start_loc:=ST_LineLocatePoint(ST_GeometryN(v_route_geom,1),v_start_pnt);
v_end_loc:=ST_LineLocatePoint(ST_GeometryN(v_route_geom,1),v_end_pnt);
-- Check if start and end points are identical to the start and end vertices
-- of the route geometry
SELECT round(ST_InterpolatePoint(ST_GeometryN(v_route_geom,1), v_start_pnt)::NUMERIC,0) INTO v_startM;
SELECT round(ST_InterpolatePoint(ST_GeometryN(v_route_geom,1), v_end_pnt)::NUMERIC,0) INTO v_endM;
-- Get the LineSubString between the two locations
v_geom_element2:= ST_LineSubstring(ST_Force2D(ST_GeometryN(v_route_geom,1)), v_start_loc, v_end_loc);
v_geom_element2:= ST_AddMeasure(v_geom_element2, 0, ST_Length(v_geom_element2));
PERFORM events.insert_split_segments(TRUE
,tmp_roads_tab
,event_name
,v_attr_val
,v_route_id
,v_route_table
,v_geom_element2
,v_route_id);
END LOOP;
ELSE
RAISE NOTICE '------------------------- ITERATION #% (%)---------------------------', cnt, event_name;
PERFORM events.create_temp_tab(tmp_roads_tab);
v_route_table:= 'temp'|| (cnt-1)::varchar || '_split_od_proskathg';
FOR event_rec IN
EXECUTE format ('SELECT id, %I::varchar, road_id, geom FROM %I WHERE road_id=$1', event_name, event_name)
USING v_route_id
LOOP
--- Assign the value of the dynamically fetched column to v_attr_val
EXECUTE format('SELECT %I FROM %I WHERE id = $1', event_name, event_name)
INTO v_attr_val
USING event_rec.id;
-- Get the route geometries from the temp (now route) table where the event intersects with
-- the route BUT, exclude the segments that intersect with only at the start or end of the event
FOR temp_route_rec IN
EXECUTE format ('SELECT n.id, from_route_id, from_route_id_src, n.geom
FROM %I n JOIN %I e ON ST_Intersects(n.geom, e.geom)
WHERE ST_Intersects(ST_EndPoint(ST_GeometryN(n.geom,1)), ST_StartPoint(e.geom))=FALSE
AND ST_Intersects(ST_StartPoint(ST_GeometryN(n.geom,1)), ST_EndPoint(e.geom))=FALSE
AND e.id=$1', v_route_table, event_name)
USING event_rec.id
LOOP
-- Get the start and points of the event
v_start_pnt:= ST_StartPoint(ST_GeometryN(event_rec.geom,1));
v_end_pnt:= ST_EndPoint(ST_GeometryN(event_rec.geom,1));
-- Check if both start and end points are along the route geometry
IF ST_DWithin(temp_route_rec.geom, v_start_pnt, 0.001) AND ST_DWithin(temp_route_rec.geom, v_end_pnt, 0.001) THEN
--They are. Now check if its at the start or the end since if it is, we need to handle it differently
-- Get the Start and End Measures
v_startM:= round(ST_InterpolatePoint(ST_GeometryN(temp_route_rec.geom,1), v_start_pnt)::NUMERIC,0);
v_endM:= round(ST_InterpolatePoint(ST_GeometryN(temp_route_rec.geom,1), v_end_pnt)::NUMERIC,0) ;
v_start_loc:=ST_LineLocatePoint(ST_GeometryN(temp_route_rec.geom,1),v_start_pnt);
v_end_loc:=ST_LineLocatePoint(ST_GeometryN(temp_route_rec.geom,1),v_end_pnt);
IF v_startM >= 0 AND v_endM <= round(ST_Length(temp_route_rec.geom)::NUMERIC,0)THEN
-- Both start and end vertices of Event lie along segment. Take the LineSubstring between the two mwasure values
v_geom_element2:= ST_LineSubstring(ST_GeometryN(temp_route_rec.geom,1), v_start_loc, v_end_loc);
v_geom_element2:= ST_AddMeasure(v_geom_element2, 0, ST_Length(v_geom_element2));
PERFORM events.insert_split_segments(FALSE
,tmp_roads_tab
,event_name
,v_attr_val
,temp_route_rec.id
,v_route_table
,v_geom_element2
,event_rec.road_id);
END IF;
ELSIF ST_DWithin(temp_route_rec.geom, v_start_pnt, 0.001)= FALSE AND ST_DWithin(temp_route_rec.geom, v_end_pnt, 0.001)= FALSE THEN
-- Both start and end vertices of Event are beyond the bounds of the segment. Take the route geometry
v_geom_element2:= temp_route_rec.geom;
PERFORM events.insert_split_segments(FALSE
,tmp_roads_tab
,event_name
,v_attr_val
,temp_route_rec.id
,v_route_table
,v_geom_element2
,event_rec.road_id);
ELSE
-- Start vertex is along the route but end vertex is beyond bounds
IF ST_DWithin(temp_route_rec.geom, v_start_pnt, 0.001)=TRUE AND ST_DWithin(temp_route_rec.geom, v_end_pnt, 0.001)= FALSE THEN
-- Start vertex of Event is along the segment but end vertex is beyond bounds of
segment Take the substring from the start location till the end of the segment
v_start_loc:=ST_LineLocatePoint(ST_GeometryN(temp_route_rec.geom,1),v_start_pnt);
-- Take the LineSubstring from start_location till the end
v_geom_element2:= ST_LineSubstring(ST_GeometryN(temp_route_rec.geom,1), v_start_loc, 1);
-- Convert to Measured Polyline
v_geom_element2:= ST_AddMeasure(v_geom_element2, 0, ST_Length(v_geom_element2));
PERFORM events.insert_split_segments(FALSE
,tmp_roads_tab
,event_name
,v_attr_val
,temp_route_rec.id
,v_route_table
,v_geom_element2
,event_rec.road_id);
ELSIF ST_DWithin(temp_route_rec.geom, v_start_pnt, 0.001)= FALSE AND ST_DWithin(temp_route_rec.geom, v_end_pnt, 0.001)= TRUE THEN
-- Start vertex of Event is beyond the bounds of the segment but end vertex is along a segment. Take the substring from the beginning till the end location
v_end_loc:=ST_LineLocatePoint(ST_GeometryN(temp_route_rec.geom,1),v_end_pnt);
-- Take the LineSubstring between the start (0) and end location
v_geom_element2:= ST_LineSubstring(ST_GeometryN(temp_route_rec.geom,1), 0, v_end_loc);
-- Convert to Measured Polyline
v_geom_element2:= ST_AddMeasure(v_geom_element2, 0, ST_Length(v_geom_element2));
PERFORM events.insert_split_segments(FALSE
,tmp_roads_tab
,event_name
,v_attr_val
,temp_route_rec.id
,v_route_table
,v_geom_element2
,event_rec.road_id);
END IF;
END IF;
END LOOP;
END LOOP;
END IF;
IF ARRAY_LENGTH(events,1)= cnt THEN
--Processed all events. Create the output table
EXECUTE format ('DROP TABLE IF EXISTS public.%I', v_output_tab);
EXECUTE FORMAT ('CREATE TABLE public.%I AS SELECT * FROM public.%I',v_output_tab, tmp_roads_tab);
END IF;
-- Increase counter
cnt:=cnt+1;
END LOOP;
END;
$function$;
Congratulations if you made it that far! If you did, and you ran the split_roads function you should see something like this (with some offsets applied on the layers). The blue line is our original network, the red, green, and orange lines are the slope15, road_width, and ortho50 events respectively and the pink line is the resulting split network. You will notice that each attribute now has the corresponding value of the event that it was split on.
That’s all for now. I hope it won’t take me another 3 years to write the next post!
Mapoholic 