-- Ensure pgrouting is installed
CREATE EXTENSION IF NOT EXISTS pgrouting;





-- Closest points on road to the STM point
SELECT r.osm_id, r.highway,
		ST_LineLocatePoint(r.way, np.way),
		ST_LineInterpolatePoint(r.way, ST_LineLocatePoint(r.way, np.way)),
		np.way <-> r.way AS distance_from_start,
		np.way, r.way
	FROM golden.natural_point np 
	INNER JOIN golden.road_line r ON ST_DWithin(np.way, r.way, 50)
	WHERE np.natural = 'peak'
		AND np.name = 'South Table Mountain'
	ORDER BY np.way <-> r.way 
;



-- Create my_route table
--   with start point for routing
DROP TABLE IF EXISTS my_route;
CREATE TABLE my_route AS 
SELECT np.way AS start_point
	FROM golden.natural_point np 
	WHERE np.name = 'South Table Mountain'
;





--  Update my_route to add end point column
ALTER TABLE my_route
	ADD end_point GEOMETRY (POINT, 3857);

-- Set end point
UPDATE my_route r
		SET end_point = ST_Centroid(b.way)
	FROM golden.building_polygon b
	WHERE b.name = 'Pangea Coffee Roasters'
;




-- Visualize start/end points
SELECT * FROM my_route;







-- Costs models for Routing:
--   Time, Distance, Scenic 
--   Use helper table for speed costs and filtering
SELECT ld.subclass, r.* 
	FROM pgosm.routable r
	INNER JOIN pgosm.layer_detail ld ON ld.code = r.code
;





/*
 * Roads need prep for routing
 */

DROP TABLE IF EXISTS pgosm.routing_roads;


-- Create roads subset filtered to sub-region.
---  Select basic attributes
CREATE TABLE pgosm.routing_roads AS
SELECT r.osm_id AS id,
		r.name, r.traffic, r.highway, r.ref, r.code,
        ST_Length(r.way)::DOUBLE PRECISION AS cost_length,
        r.way AS the_geom
    FROM ( -- Create single Bounding Box for a specific region
    	SELECT ST_Envelope(ST_Union(way)) AS way
            FROM golden.boundary_polygon bp 
            WHERE bp.name = 'Golden'
        ) b
    INNER JOIN golden.road_line r 
    	ON r.way && b.way  
    INNER JOIN pgosm.routable rbl 
    	ON r.code = rbl.code
    		 -- Routes for non-motorized traffic
    		AND rbl.route_foot
    		AND NOT rbl.route_motor
;



CREATE INDEX GIX_pgosm_routing_roads ON pgosm.routing_roads
    USING GIST (the_geom);

   
   

-- Ensure pgrouting-generated tables don't already exist
DROP TABLE IF EXISTS pgosm.routing_roads_noded;
DROP TABLE IF EXISTS pgosm.routing_roads_noded_vertices_pgr;





-- Run pgrouting steps on custom roads layer
--    Creates <table>_noded
--    Expensive step!
SELECT pgr_nodeNetwork('pgosm.routing_roads', 1);



-- Creates <table>_noded_vertices_pgr
--    Expensive step!
SELECT pgr_createTopology('pgosm.routing_roads_noded', 0.1);


-- Analyzes <table>_noded_vertices_pgr
--   Populates `cnt` and `chk`
--   chk > 0 indicates vertix might have problem
SELECT pgr_analyzegraph('pgosm.routing_roads_noded', 0.1);





-- Noded versions have more rows
SELECT 'original' AS src, COUNT(*) FROM pgosm.routing_roads
UNION
SELECT 'noded' AS src, COUNT(*) FROM pgosm.routing_roads_noded
;



-- old_id links to osm_id
SELECT *
	FROM pgosm.routing_roads_noded
	WHERE old_id = 29276956;


SELECT osm_id, highway, way
	FROM golden.road_line 
	WHERE osm_id = 29276956;


-- Update <table>_noded with extra attributes for routing
ALTER TABLE pgosm.routing_roads_noded ADD cost_length DOUBLE PRECISION NULL;
ALTER TABLE pgosm.routing_roads_noded ADD cost_minutes DOUBLE PRECISION NULL;
ALTER TABLE pgosm.routing_roads_noded ADD code TEXT NULL;






-- Update noded with costs and sub-class code
UPDATE pgosm.routing_roads_noded rn
    SET code = r.code,
    	cost_length = ST_Length(rn.the_geom)::DOUBLE PRECISION,
        -- Math to calculate cost (minutes) based on length (meters)
        --   and speed (MPH) for given road classification
        cost_minutes = ST_Length(rn.the_geom)::DOUBLE PRECISION
			/ (rbl.max_speed / 60 * 1.609344 * 1000)
    FROM pgosm.routing_roads r
    INNER JOIN pgosm.routable rbl ON r.code = rbl.code
    -- Maps to original OSM data osm_id
    WHERE rn.old_id = r.id
;






-- Clear out unroutable roads
--   Check your data before deleting it! ;-)
DELETE FROM pgosm.routing_roads_noded WHERE source IS NULL;
DELETE FROM pgosm.routing_roads_noded WHERE cost_minutes IS NULL;





/*
 * Roads are now routing-ready!
 */






-- Visualize options for Route End Point
--    <table>_noded_vertices_pgr has the points our
--    road data can route to
SELECT r.id, a.end_point, r.the_geom
	FROM my_route a
	INNER JOIN pgosm.routing_roads_noded_vertices_pgr r
		ON ST_DWithin(a.end_point, r.the_geom, 200)
	ORDER BY a.end_point <-> r.the_geom ASC
	LIMIT 15
;





-- Add columns for start/end nodes (noded vertices for roads)
ALTER TABLE my_route ADD start_node_id INT;
ALTER TABLE my_route ADD end_node_id INT;



SELECT * FROM my_route;







-- Update my_route with Start Node ID
UPDATE my_route 
	SET start_node_id = o.id
	FROM (
		SELECT r.id 
			FROM my_route a
			INNER JOIN pgosm.routing_roads_noded_vertices_pgr r
				ON ST_DWithin(a.start_point, r.the_geom, 200)
			ORDER BY a.start_point <-> r.the_geom ASC
			LIMIT 1
		) o
;
		
		
-- Update my_route with End Node ID
UPDATE my_route 
	SET end_node_id = o.id
	FROM (
		SELECT r.id 
			FROM my_route a
			INNER JOIN pgosm.routing_roads_noded_vertices_pgr r
				ON ST_DWithin(a.end_point, r.the_geom, 200)
			ORDER BY a.end_point <-> r.the_geom ASC
			LIMIT 1
		) o
;




-- Join my_route to routing start/end details
SELECT a.*, b.the_geom AS route_start, c.the_geom AS route_end
	FROM my_route a
	INNER JOIN pgosm.routing_roads_noded_vertices_pgr b 
		ON a.start_node_id = b.id
	INNER JOIN pgosm.routing_roads_noded_vertices_pgr c
		ON a.end_node_id = c.id
	;




-- pgr_dijkstra() to route
--   Provides step-by-step route with costs
SELECT p.start_point, p.end_point, p.start_node_id, p.end_node_id,
		ld.subclass, rbl.max_speed, c.cost, c.agg_cost, n.the_geom
    FROM my_route p
    INNER JOIN pgr_dijkstra('SELECT id, source, target,
                    cost_minutes AS cost
                FROM pgosm.routing_roads_noded ',
            p.start_node_id,
            p.end_node_id,
            False) c
        ON True
    INNER JOIN pgosm.routing_roads_noded n ON c.edge = n.id
    INNER JOIN pgosm.layer_detail ld ON n.code = ld.code
    INNER JOIN pgosm.routable rbl ON ld.code = rbl.code
	ORDER BY path_seq
;



ALTER TABLE my_route ADD route GEOMETRY(MULTILINESTRING, 3857);
ALTER TABLE my_route ADD cost_minutes NUMERIC(12,2);
ALTER TABLE my_route ADD cost_length NUMERIC(12,2);





-- Wrap into CTE for UPDATE
WITH a AS (
SELECT SUM(c.cost) AS agg_cost, 
		ST_Collect(n.the_geom) AS way
    FROM my_route p
    INNER JOIN pgr_dijkstra('SELECT id, source, target,
                    cost_minutes AS cost
                FROM pgosm.routing_roads_noded ',
            p.start_node_id,
            p.end_node_id,
            False) c
        ON True
    INNER JOIN pgosm.routing_roads_noded n ON c.edge = n.id
)
UPDATE my_route 
	SET cost_minutes = agg_cost,
		route = way,
		cost_length = ST_Length(way)
	FROM a
;



SELECT * 
	FROM my_route;







SELECT * FROM public.layer_styles;
-- Visualize in QGIS