road_network_extraction_with_pr_routing_and_plpython

--------------------------
-- Rémi Cura, Thales IGN, 2015
-- this script uses pgrouting to leverage connectivity of patches
-- useful for the journal article about point cloud server
--
-- NOTE: depends on rc_lib extension  : https://github.com/Remi-C/PPPP_utilities
---------------------------


--create a schema, set search_path
/*
CREATE SCHEMA IF NOT EXISTS patch_connectivity ; 
SET search_path to patch_connectivity,patch_connectivity_clustering,  benchmark_cassette_2013,tmob_20140616, test_grouping, public , topology;

--creating the table to hold nodes
DROP TABLE IF EXISTS patch_node_table CASCADE ;
CREATE TABLE patch_node_table (
 node_id SERIAL PRIMARY KEY
 , patch_centroid geometry(pointZ, 932011 ) 
 , patch_height float
 , spatial_size float
) ; 
TRUNCATE patch_node_table CASCADE ;
CREATE INDEX ON patch_node_table USING GIST(patch_centroid) ; 
CREATE INDEX ON patch_node_table USING GIST(patch_centroid gist_geometry_ops_nd) ; 
CREATE INDEX ON patch_node_table (spatial_size ) ; 
--filling the table with all patches
TRUNCATE patch_node_table CASCADE ; 
INSERT INTO patch_node_table 
SELECT DISTINCT ON (floor(ST_X(center)  ), floor(ST_Y(center)    ),floor( (upper(avg_z) + lower(avg_z)) / 2.0 ) ) 
	gid,
	ST_SetSRID(ST_MakePoint(st_x(center), st_y(center), pc_patchavg(patch,'z')),932011)
	-- ST_SetSRID(ST_MakePoint(round(ST_X(center)), round(ST_Y(center)),round(height_above_laser ) )  , 932011 ) 
	,floor( (upper(avg_z) + lower(avg_z)) / 2.0 ) 
FROM riegl_pcpatch_space_int_proxy  AS prox LEFT OUTER  JOIN riegl_pcpatch_space_int USING (gid), ST_Centroid(geom) AS center
	, -- ST_SetSRID(ST_MAkePoint(1907.77,21141.02),932011 ) as ref_point
	ST_SetSRID(ST_MAkePoint(650907.6841,6861141.0047),931008 ) as ref_point
	
WHERE -- ST_DWithin(geom, ref_point,30)  
	--AND
	 num_points > 100
	 AND ST_Area(prox.geom) > 0.6
	--	AND gid %4 = 0
ORDER BY floor(ST_X(center)  ), floor(ST_Y(center)    ),floor( (upper(avg_z) + lower(avg_z)) / 2.0 ),num_points DESC


TRUNCATE patch_node_table CASCADE; 
ALTER SEQUENCE patch_node_table_node_id_seq RESTART WITH 1;
INSERT INTO patch_node_table
SELECT   DISTINCT ON (floor(ST_X(center)  ), floor(ST_Y(center)    ),floor( (upper(avg_z) + lower(avg_z)) / 2.0 ) ) 
	row_number()  over() -- gid 
	, ST_SetSRID(ST_MakePoint(st_x(center), st_y(center), (upper(avg_z) + lower(avg_z)) / 2.0 ),932011)
	, (upper(avg_z) + lower(avg_z)) / 2.0 
	FROM ST_GeomFromText('POINT(650967.0 6861577.5)',931008) as ref_point, 
		tmob_20140616.riegl_pcpatch_space_int_proxy, ST_Centroid(geom) AS center
WHERE  ST_DWithin(geom, ref_point,50)  AND
	 num_points > 100

		-- AND --(ST_Area(geom) > 0.95 OR (upper(avg_z) - lower(avg_z)) < 0.3) 
		--	  upper(avg_z) - lower(avg_z)< 0.4
			-- (upper(avg_z) - lower(avg_z)) * (1-ST_Area(geom)) < 0.1
		-- AND gid % 2 = 0
ORDER  BY floor(ST_X(center)  ), floor(ST_y(center)    ), floor( (upper(avg_z) + lower(avg_z)) / 2.0) , num_points ASC
 	 LIMIT 20000

	 SELECT *
	 FROM tmob_20140616.riegl_pcpatch_space_int_proxy
	 LIMIT 1 
-- WHERE gid < 1000; 
-- UPDATE patch_node_table SET spatial_size = 1 
 

--filling with patch of varying size
TRUNCATE patch_node_table  CASCADE; 
INSERT INTO patch_node_table 
SELECT gid ,-- ST_SetSRID(ST_MakePoint(ST_x(center),  ST_y(center), avg_z),932011) 
		ST_SetSRID(
		ST_MakePoint(
			  ST_X(center)   
			,  ST_Y(center)    
			,  avg_z   
			)
			, 932011 ) 
	 --, ST_Force3D(ST_SnaptoGrid(center, spatial_size ) ) 
	, NULL AS height_above_laser
	, spatial_size
FROM copy_bench, ST_Centroid(patch::geometry) AS center

SELECT *
FROM copy_bench, ST_Centroid(patch::geometry) AS center
WHERE avg_z IS NULL OR ST_IsEmpty(center)  

SELECT count(*)  
FROM patch_node_table
LIMIT 100 ; 

--computing an edge table 
DROP TABLE IF EXISTS edge_table ; 
CREATE TABLE edge_table  (
	edge_id SERIAL PRIMARY KEY
	, source INT REFERENCES patch_node_table (node_id)
	, target INT REFERENCES patch_node_table (node_id)
	, tcost float
	, geom geometry(linestringZ, 932011) 
) ; 
CREATE INDEX ON edge_table USING GIST(geom gist_geometry_ops_nd)  ; 
CREATE INDEX ON edge_table USING GIST(geom )  ; 

TRUNCATE edge_table ;

SELECT count(*)
FROM edge_table
WHERE tcost = 0 ;  


INSERT INTO edge_table (source,target, tcost, geom)
SELECT et.target, et.source, et.tcost, ST_Reverse(geom)
FROM edge_table as et

SELECT *
FROM edge_table 
LIMIT 1

-- for ign server
TRUNCATE edge_table ; 
WITH node_connected AS ( --721466
	SELECT  p1.node_id AS node_1, p2.node_id AS node_2 
		, ST_3DDistance(p1.patch_centroid, p2.patch_centroid) AS tcost
			--CASE WHEN po1.dominant_simplified_class != po2.dominant_simplified_class THEN 100 ELSE ST_3DDistance(p1.patch_centroid, p2.patch_centroid) END AS tcost
		, ST_SetSRID(ST_MakeLine(p1.patch_centroid, p2.patch_centroid),932011) as line
	FROM   patch_node_table AS p1 
		, patch_node_table AS p2 
	WHERE p1.node_id < p2.node_id
		AND ST_3DDWithin(p1.patch_centroid, p2.patch_centroid,1.5) = TRUE -- SELECT sqrt(3)+0.1 --SELECT sqrt(0.5+0.5^2) -- SELECT sqrt(2) 
		--AND ST_DWithin(p1.patch_centroid, ref_point,120) 
		--AND ST_DWithin(p2.patch_centroid,ref_point,120) 
		 -- 1.1 : 4 connectivity
		 -- 1.5 : 8 connectivity
		 -- 1.9 : 3D 8 connectivity
		  --AND p1.node_id < 100
)
INSERT INTO edge_table (source, target, tcost,geom) SELECT node_1, node_2, tcost , line
FROM node_connected ; 

--removing nodes that are in no edges
DELETE FROM patch_node_table AS pn WHERE NOT EXISTS (SELECT 1 FROM edge_table AS et WHERE pn.node_id = et.source OR pn.node_id = et.target ) ; 


--for varying size
TRUNCATE edge_table ; 
WITH node_connected AS ( --721466
	SELECT  p1.node_id AS node_1, p2.node_id AS node_2 
		, ST_3DDistance(p1.patch_centroid, p2.patch_centroid) AS tcost
			--CASE WHEN po1.dominant_simplified_class != po2.dominant_simplified_class THEN 100 ELSE ST_3DDistance(p1.patch_centroid, p2.patch_centroid) END AS tcost
		, ST_SetSRID(ST_MakeLine(p1.patch_centroid, p2.patch_centroid),932011) as line
	FROM ST_SetSRID(ST_MAkePoint(1907.77,21141.02),932011 ) as ref_point,  patch_node_table AS p1
		LEFT OUTER JOIN riegl_pcpatch_proxy AS po1 ON (po1.gid = p1.node_id)
		, patch_node_table AS p2
		LEFT OUTER JOIN riegl_pcpatch_proxy AS po2 ON (po2.gid = p2.node_id)
	WHERE p1.node_id < p2.node_id
		AND ST_3DDWithin(p1.patch_centroid, p2.patch_centroid,2  ) = TRUE -- SELECT sqrt(3)+0.1 --SELECT sqrt(0.5+0.5^2) -- SELECT sqrt(2)
		AND ST_3DDWithin(p1.patch_centroid, p2.patch_centroid,sqrt( (p1.spatial_size + p2.spatial_size + p1.spatial_size * p2.spatial_size)) ) = TRUE 
		AND ST_DWithin(p1.patch_centroid, ref_point,30) 
		AND ST_DWithin(p2.patch_centroid,ref_point,30) 
		 -- 1.1 : 4 connectivity
		 -- 1.5 : 8 connectivity
		 -- 1.9 : 3D 8 connectivity
)
INSERT INTO edge_table (source, target, tcost,geom) SELECT node_1, node_2, tcost , line
FROM node_connected ; 

-- for local benchmark
TRUNCATE edge_table ; 
WITH node_connected AS ( --721466
	SELECT  p1.node_id AS node_1, p2.node_id AS node_2 
		, ST_3DDistance(p1.patch_centroid, p2.patch_centroid) AS tcost
			--CASE WHEN po1.dominant_simplified_class != po2.dominant_simplified_class THEN 100 ELSE ST_3DDistance(p1.patch_centroid, p2.patch_centroid) END AS tcost
		, ST_SetSRID(ST_MakeLine(p1.patch_centroid, p2.patch_centroid),932011) as line
	FROM ST_SetSRID(ST_MAkePoint(1907.77,21141.02),932011 ) as ref_point, patch_node_table AS p1
		LEFT OUTER JOIN copy_bench AS po1 ON (po1.gid = p1.node_id)
		, patch_node_table AS p2
		LEFT OUTER JOIN copy_bench AS po2 ON (po2.gid = p2.node_id)
	WHERE p1.node_id < p2.node_id
		AND ST_3DDWithin(p1.patch_centroid, p2.patch_centroid,3 ) = TRUE -- SELECT sqrt(3)+0.1 --SELECT sqrt(0.5+0.5^2) -- SELECT sqrt(2)
		-- AND ST_3DDWithin(p1.patch_centroid, p2.patch_centroid,sqrt(sqrt( p1.spatial_size* p2.spatial_size) + p1.spatial_size*p2.spatial_size))  = TRUE 
		--AND ST_DWithin(p1.patch_centroid, ref_point,20) 
		--	 AND ST_DWithin(p2.patch_centroid, ref_point,20) 
		 -- 1.1 : 4 connectivity
		 -- 1.5 : 8 connectivity
		 -- 1.9 : 3D 8 connectivity
)
INSERT INTO edge_table (source, target, tcost,geom) SELECT node_1, node_2, tcost , line
FROM node_connected ; 

SELECT ST_SRID(patch_centroid)
FROM patch_node_table
LIMIT 1 

-----------------
-- exporting node and edge for cloud compare
------------------
 pgsql2shp -f patch_node_table -h localhost -p 5433 -u postgres -P postgres -g patch_centroid -r test_pointcloud patch_connectivity.patch_node_table
 pgsql2shp -f patch_node_table -h 172.16.3.50 -p 5432 -u postgres -P postgres -g patch_centroid -r test_pointcloud patch_connectivity.patch_node_table

 pgsql2shp -f edge_table -h 172.16.3.50 -p 5432 -u postgres -P postgres -g geom -r test_pointcloud patch_connectivity.edge_table

 COPY (
	SELECT  ST_X(point) AS x , ST_Y(point) AS Y, ST_Z(point) AS Z  
	FROM ST_GeomFromText('POINT(650967.0 6861577.5)',931008) as ref_point, 
		tmob_20140616.riegl_pcpatch_space_int_proxy as prox
			LEFT OUTER JOIN tmob_20140616.riegl_pcpatch_space_int as pa USING (gid)
			, rc_exploden_grid(patch,1000,0.025) as pt
			, CAST(pt AS geometry) AS point
		WHERE  ST_DWithin(prox.geom, ref_point,50) 
 
) TO '/tmp/export_points.csv' WITH CSV HEADER ;
------------------

-- creating pgrouting extension
CREATE EXTENSION IF NOT EXISTS pgrouting  ; 

--finding 2 patches (start/end)
DROP TABLE IF EXISTS k_shortest_path ;
CREATE TABLE k_shortest_path AS  
	WITh start_node AS (
		SELECT node_id as snode
		FROM patch_node_table
			--, ST_SetSRID(ST_MakePoint(1900.65,21124.27,45),932011) AS start_point
			 --,ST_SetSRID(ST_MakePoint( 1899.89,21156.33,45),932011) AS start_point
			 --
			 ,ST_SetSRID(ST_MakePoint( 650900.649778,6861124.270031,45),932011) AS start_point
		WHERE ST_3DDWithin(patch_centroid ,start_point , 2 ) 
			ORDER BY ST_3DDistance(patch_centroid ,start_point ) ASC
			LIMIT 1 
	)
	, end_node AS (
		SELECT node_id AS enode
		FROM patch_node_table 
			-- , ST_SetSRID(ST_MakePoint(1910.13,21156.33,45),932011) AS end_point
			, ST_SetSRID(ST_MakePoint(650910.129758,6861156.329822,45),932011) AS end_point
		WHERE ST_3DDWithin(patch_centroid ,end_point , 2 ) 
			ORDER BY ST_3DDistance(patch_centroid ,end_point ) ASC
			LIMIT 1
	)
	  , pgrouting_results AS (
		 SELECT   f.* --seq,  path_id, path_seq,node, edge AS edge_id, agg_cost, cost
		 FROM start_node, end_node 
			--, pgr_ksp(
			--   'SELECT edge_id AS id, source, target, tcost AS cost, tcost AS reverse_cost FROM edge_table',
			--    snode,enode
			--    , 1
			--    ,directed:=false
			-- ) 
		 ,pgr_astar(
			   'SELECT edge_id AS id, source, target, tcost AS cost,St_x(spoint) AS x1, St_x(epoint) AS x2, St_y(spoint)AS y1, St_y(epoint) AS y2
			    FROM edge_table, ST_StartPoint(geom) as spoint, ST_EndPoint(geom) as epoint ',
			    snode,enode
			    ,directed:=false
			    ,has_reverse_cost:=false  
			 ) as f
	 )
	 SELECT seq AS uid , path_id, path_seq, edge_id, agg_cost, geom::geometry(linestringZ,932011) AS geom
		, ST_Force2D(geom)::geometry(linestring,932011) AS geom2D
	 FROM pgrouting_results AS pg
		NATURAL JOIN edge_table AS et1   ; 
	CREATE INDEX ON k_shortest_path  (uid) ; 
CREATE INDEX ON k_shortest_path USING GIST(geom) ; 
CREATE INDEX ON k_shortest_path USING GIST(geom2D) ; 


--------using astar------------

--finding 2 patches (start/end)
DROP TABLE IF EXISTS astar;
CREATE TABLE astar AS  
	WITh start_node AS (
		SELECT node_id as snode
		FROM patch_node_table 
			 ,ST_SetSRID(ST_MakePoint( 650900.649778,6861124.270031,45),932011) AS start_point
		WHERE ST_3DDWithin(patch_centroid ,start_point , 2 ) 
			ORDER BY ST_3DDistance(patch_centroid ,start_point ) ASC
			LIMIT 1 
	)
	, end_node AS (
		SELECT node_id AS enode
		FROM patch_node_table  
			, ST_SetSRID(ST_MakePoint(650910.129758,6861156.329822,45),932011) AS end_point
		WHERE ST_3DDWithin(patch_centroid ,end_point , 2 ) 
			ORDER BY ST_3DDistance(patch_centroid ,end_point ) ASC
			LIMIT 1
	)
	 , pgrouting_results AS (
		 SELECT   seq, id1 AS node_id ,id2 AS edge_id ,cost -- ,  path_id, path_seq,node, edge AS edge_id, agg_cost, cost
		 FROM start_node, end_node  
			 ,pgr_astar(
				   'SELECT edge_id AS id, source, target, tcost AS cost,St_x(spoint) AS x1,  St_y(spoint)AS y1, St_x(epoint) AS x2, St_y(epoint) AS y2
				    FROM edge_table, ST_StartPoint(geom) as spoint, ST_EndPoint(geom) as epoint ',
				    snode,enode
				    ,directed:=false
				    ,has_reverse_cost:=false  
				 ) as f
	 )
	 SELECT seq AS uid , node_id, edge_id, cost, ST_SetSRID(geom,931008)::geometry(linestringZ,931008) AS geom
		, ST_SetSRID(ST_Force2D(geom),931008) ::geometry(linestring,931008) AS geom2D
	 FROM pgrouting_results AS pg
		NATURAL JOIN edge_table AS et1   ; 
	CREATE INDEX ON astar  (uid) ; 
CREATE INDEX ON astar USING GIST(geom) ; 
CREATE INDEX ON astar USING GIST(geom2D) ; 



DROP TABLE IF EXISTS astar_export ;
CREATE TABLE astar_export AS
SELECT  1::int as id,  ST_Translate(ST_LineMerge(ST_Collect(geom ORDER BY uid ASC)) ,0,0,0 )::geometry(LinestringZ,931008) AS geom
FROM astar ; 


-- pgsql2shp -f shortest_path -h 172.16.3.50 -p 5432 -u postgres -P postgres -g geom -r test_pointcloud patch_connectivity.astar_export


--exporting point cloud to visualize
COPY (
	 SELECT st_x(pt) AS x, st_y(pt) AS y, st_z(pt) AS z
	 FROM patch_node_table LEFT OUTER JOIN tmob_20140616.riegl_pcpatch_space_int  as pa ON (node_id = gid)
		, pc_explode(patch) AS point, CAST(point AS geometry) as pt 
) TO '/ExportPointCloud/testing_distance_variation.csv' WITH CSV HEADER ;

--export point cloud that is near the astar path 
COPY (
	 SELECT st_x(pt) AS x, st_y(pt) AS y, st_z(pt) AS z
	 FROM astar_export AS ast, patch_node_table AS pn LEFT OUTER JOIN tmob_20140616.riegl_pcpatch_space_int  as pa ON (node_id = gid)
		, rc_exploden_grid(patch,1000, 0.04) AS point, CAST(point AS geometry) as pt 
	WHERE ST_3DDWithin(ST_SetSRID(ast.geom,932011),pn.patch_centroid,2) 
		AND ST_3DDWithin( ast.geom ,pt,2) 
) TO '/ExportPointCloud/testing_distance_variation_only_point_on_path.csv' WITH CSV HEADER ;


--now, again using pgrouting to compute the shortest path with points 
--creating the table to hold nodes whch are points
DROP TABLE IF EXISTS point_node_table CASCADE ;
CREATE TABLE point_node_table (
 node_id SERIAL PRIMARY KEY
 , patch_centroid geometry(pointZ, 931008 ) 
 , patch_height float
 , spatial_size float
) ;  
CREATE INDEX ON point_node_table USING GIST(patch_centroid) ; CREATE INDEX ON point_node_table USING GIST(patch_centroid gist_geometry_ops_nd) ; CREATE INDEX ON point_node_table (spatial_size ) ; 

TRUNCATE point_node_table CASCADE ;
INSERT INTO point_node_table 
	SELECT row_number() over() AS node_id
		, pt as patch_centroid
		, ST_Z(pt) AS patch_height
		, 1 AS spatila_size
	FROM astar_export AS ast, patch_node_table AS pn LEFT OUTER JOIN tmob_20140616.riegl_pcpatch_space_int  as pa ON (node_id = gid)
		, rc_exploden_grid(patch,1000, 0.04) AS point, CAST(point AS geometry) as pt 
	WHERE ST_3DDWithin(ST_SetSRID(ast.geom,932011),pn.patch_centroid,2) 
		AND ST_3DDWithin( ast.geom ,pt,2) ; 
		
DROP TABLE IF EXISTS point_edge_table ; 
CREATE TABLE point_edge_table  (
	edge_id SERIAL PRIMARY KEY
	, source INT REFERENCES point_node_table (node_id)
	, target INT REFERENCES point_node_table (node_id)
	, tcost float
	, geom geometry(linestringZ, 932011) 
) ; 
CREATE INDEX ON point_edge_table USING GIST(geom gist_geometry_ops_nd)  ; 
CREATE INDEX ON point_edge_table USING GIST(geom )  ; 

TRUNCATE point_edge_table ; 
WITH node_connected AS ( --721466
	SELECT  p1.node_id AS node_1, p2.node_id AS node_2 
		, ST_3DDistance(p1.patch_centroid, p2.patch_centroid) AS tcost 
		, ST_SetSRID(ST_MakeLine(p1.patch_centroid, p2.patch_centroid),932011) as line
	FROM   point_node_table AS p1 
		, point_node_table AS p2 
	WHERE p1.node_id < p2.node_id
		AND ST_3DDWithin(p1.patch_centroid, p2.patch_centroid,0.1) = TRUE 
)
INSERT INTO point_edge_table (source, target, tcost,geom) SELECT node_1, node_2, tcost , line
FROM node_connected ; 

-- now getting the shortest path using points :

DROP TABLE IF EXISTS points_astar;
CREATE TABLE astar AS  
	WITh start_node AS (
		SELECT node_id as snode
		FROM point_node_table 
			 ,ST_SetSRID(ST_MakePoint( 650900.649778,6861124.270031,45),931008) AS start_point
		WHERE ST_3DDWithin(patch_centroid ,start_point , 2 ) 
			ORDER BY ST_3DDistance(patch_centroid ,start_point ) ASC
			LIMIT 1 
	)
	, end_node AS (
		SELECT node_id AS enode
		FROM point_node_table  
			, ST_SetSRID(ST_MakePoint(650910.129758,6861156.329822,45),931008) AS end_point
		WHERE ST_3DDWithin(patch_centroid ,end_point , 2 ) 
			ORDER BY ST_3DDistance(patch_centroid ,end_point ) ASC
			LIMIT 1
	)
	-- , pgrouting_results AS (
		 SELECT   seq, id1 AS node_id ,id2 AS edge_id ,cost -- ,  path_id, path_seq,node, edge AS edge_id, agg_cost, cost
		 FROM start_node, end_node  
			 ,pgr_astar(
				   'SELECT edge_id AS id, source, target, tcost AS cost,St_x(spoint) AS x1,  St_y(spoint)AS y1, St_x(epoint) AS x2, St_y(epoint) AS y2
				    FROM point_edge_table, ST_StartPoint(geom) as spoint, ST_EndPoint(geom) as epoint ',
				    snode,enode
				    ,directed:=false
				    ,has_reverse_cost:=false  
				 ) as f
	 )
	 SELECT seq AS uid , node_id, edge_id, cost, ST_SetSRID(geom,931008)::geometry(linestringZ,931008) AS geom
		, ST_SetSRID(ST_Force2D(geom),931008) ::geometry(linestring,931008) AS geom2D
	 FROM pgrouting_results AS pg
		NATURAL JOIN point_edge_table AS et1   ; 
	CREATE INDEX ON point_astar  (uid) ; 
CREATE INDEX ON point_astar USING GIST(geom) ; 
CREATE INDEX ON point_astar USING GIST(geom2D) ; 



DROP TABLE IF EXISTS astar_export ;
CREATE TABLE astar_export AS
SELECT  1::int as id,  ST_Translate(ST_LineMerge(ST_Collect(geom ORDER BY uid ASC)) ,0,0,0 )::geometry(LinestringZ,931008) AS geom
FROM astar ; 

--------------------------------


--testing a query with 2 example nodes :  756 1089
DROP TABLE IF EXISTS k_shortest_path ;
CREATE TABLE k_shortest_path AS  
WITH pgrouting_results AS (
	 SELECT   seq,  path_id, path_seq,node, edge AS edge_id, agg_cost, cost
	 FROM pgr_ksp(
		   'SELECT edge_id AS id, source, target, tcost AS cost, tcost AS reverse_cost FROM edge_table',
		   -- 753, 4005--north to south, ground
		  --  15008, 20884 --wall to wall
		    10631,17339
		    , 10 
		    ,directed:=false
		 )
 )
 SELECT seq AS uid , path_id, path_seq, edge_id, agg_cost, geom::geometry(linestringZ,932011) AS geom
 FROM pgrouting_results AS pg
	NATURAL JOIN edge_table AS et1   ; 
CREATE INDEX ON k_shortest_path USING GIST(geom) ; 




--exporting point cloud
-- /media/sf_E_RemiCura/PROJETS/articles_ISPRS_Geospatial_Week/experiments/pgrouting$ pgsql2shp -f shortest_path -h localhost -p 5433 -u postgres -P postgres -g geom -r test_pointcloud patch_connectivity.k_shortest_path_export

DROP TABLE IF EXISTS k_shortest_path_export ;
CREATE TABLE k_shortest_path_export AS
SELECT path_id ,  ST_Translate(ST_LineMerge(ST_Collect(geom ORDER BY uid ASC)) ,0,0,0 )::geometry(LinestringZ,932011) AS geom
FROM k_shortest_path
GROUP BY path_id ; 




--computing the isochrone
	--1902.98,21263.10
DROP TABLE IF EXISTS isochrone ; 
CREATE TABLE  isochrone  AS
WITH isochrone_result AS (
	SELECT * -- seq, id1 AS node, id2 AS edge, cost, geom
	  FROM pgr_drivingdistance( 
	    'SELECT edge_id AS id, source, target, tcost AS cost 
	    FROM ST_SetSRID(ST_MakePoint(1907.03,21197.88,45),932011)as tpoint, edge_table  
	    WHERE ST_3DDWithin(tpoint, geom, 50) = TRUE
		 ',
	    15331, 50 , directed := false
	  ) as di  
)
SELECT seq, node_id, edge, agg_cost, patch_centroid::geometry(pointZ, 932011)
FROM isochrone_result AS iso
	LEFT OUTER JOIN patch_node_table AS pn ON (pn.node_id = iso.node) ;


DROP TABLE IF EXISTS isochrone_export  ;
CREATE TABLE  isochrone_export  AS 
SELECT node_id, agg_cost, patch_centroid::geometry(pointZ, 932011)
FROM isochrone  ; 

--exporting the points as csv for cloudcompare
COPY (
	SELECT DISTINCT ON (node_id ) X,Y,Z,agg_cost, dominant_simplified_class
	FROM (
	SELECT node_id,  ST_X(patch_centroid) AS X, ST_Y(patch_centroid) AS Y, ST_Z(patch_centroid) +40.5 AS Z, agg_cost , dominant_simplified_class
	FROM isochrone_export AS pt
		LEFT OUTER JOIN riegl_pcpatch_proxy AS rpp ON (rpp.gid = pt.node_id)
	UNION ALL 
	SELECT node_id,  ST_X(patch_centroid) AS X, ST_Y(patch_centroid) AS Y, ST_Z(patch_centroid) +40.5 AS Z, -1 AS agg_cost ,dominant_simplified_class
	FROM patch_node_table  AS pt
		LEFT OUTER JOIN riegl_pcpatch_proxy AS rpp ON (rpp.gid = pt.node_id)
	WHERE ST_3DDWithin(ST_SetSRID(ST_MakePoint(1907.03,21197.88,45),932011), patch_centroid , 50) = TRUE
		AND dominant_simplified_class  = 4 --building
	) AS sub
	ORDER BY node_id, agg_cost DESC

) TO '/media/sf_E_RemiCura/PROJETS/articles_ISPRS_Geospatial_Week/experiments/pgrouting/test_isochrone.csv' WITH CSV HEADER ;


CREATE EXTENSION IF NOT EXISTS pgrouting
-- DROP TABLE isochrone_geomdist  ;
-- CREATE TABLE isochrone_geomdist  AS 

DROP TABLE IF EXISTS all_shortest_path_3  ;
CREATE TABLE all_shortest_path_3 AS
-- all to all topological distance
WITH all_shortest_path AS (
	SELECT seq, id1, id2, cost
	FROM  pgr_apspJohnson(
	' SELECT *
	FROM (SELECT edge_id AS id, source, target, tcost AS cost  
		    FROM  edge_table UNION  ALL
	SELECT 1000000+edge_id AS id, target, source, tcost AS cost
	FROM edge_table) AS sub
	LIMIT 100'::text)-- as f , false, false ) AS f

	    )
SELECT  seq, St_SetSRID(ST_MakeLine(pn1.patch_centroid,pn2.patch_centroid),932011)::geometry(linestringZ,932011) AS geom
 FROM all_shortest_path AS pg
	LEFT OUTER JOIN patch_node_table AS pn1 ON (pn1.node_id = id1)
	LEFT OUTER JOIN patch_node_table AS pn2 ON (pn2.node_id = id2)


-- writing a function to compute shortest path cost between a node and k other nodes

--generating all pairs within 50 m of distance

DROP TABLE IF EXISTS shortest_path_DWithin ; 
CREATE TABLE shortest_path_DWithin AS 
WITH all_pairs AS (
	SELECT DISTINCT et1.node_id AS node_1 , et2.node_id AS node_2
	FROM patch_node_table AS et1, patch_node_table AS et2
	WHERE ST_DWithin(et1.patch_centroid, et2.patch_centroid, 50)  = TRUE
		AND et1.node_id < et2.node_id --7135
		AND --checking that the nodes are used in edges, else isolated nodes could be a problem
			EXISTS (SELECT 1 FROM edge_table AS et 
				WHERE (et.source = et1.node_id OR et.source = et2.node_id) ) 
		AND EXISTS (SELECT 1 FROM edge_table AS et 
				WHERE (et.target = et1.node_id OR et.target = et2.node_id) ) 
 
	ORDER BY node_1, node_2
)
 , pair_agg AS (
	SELECT node_1, array_agg(node_2 ORDER BY node_2 ASC) AS node_2s
	FROM all_pairs
	GROUP BY node_1
)
, pair_cost AS ( 
	SELECT f.*
	FROM pair_agg , pgr_kdijkstraCost(
		'SELECT edge_id AS id, source, target, tcost AS cost FROM edge_table'
		, source_vid:= node_1
		, target_vid:=node_2s
		, directed:=FALSE
		,has_reverse_cost:=FALSE) AS f
	
)
SELECT seq, id1 AS node_id1, id2 AS node_id2, cost AS tcost, St_SetSRID(ST_MakeLine(pn1.patch_centroid,pn2.patch_centroid),932011)::geometry(linestringZ,932011) AS geom
FROM pair_cost  AS pg
	LEFT OUTER JOIN patch_node_table AS pn1 ON (pn1.node_id = id1)
	LEFT OUTER JOIN patch_node_table AS pn2 ON (pn2.node_id = id2) 
WHERE cost >0; 

 
SELECT *
FROM shortest_path_DWithin
WHERE tcost < 10
LIMIT 1
 */
/*
	DROP TABLE IF EXISTS network_clustering ; 
	CREATE TABLE network_clustering AS 
		WITH preparing_input AS (
			SELECT array_agg(node_id1 ORDER BY seq) as node_id1s, array_agg(node_id2 ORDER BY seq) as node_id2s, array_agg(tcost ORDER BY seq) as tcosts 
			FROM shortest_path_DWithin
			-- WHERE tcost < 3
		)
		 , result_clustering AS (
		SELECT f.*
		FROM preparing_input,  rc_graph_clustering(node_id1s,node_id2s,tcosts,2) AS f
		)
		 , artificial_seq  AS (
			SELECT DISTINCT row_number() over(order by node_id) -1 as seq,   node_id 
			FROM (SELECT node_id1 AS node_id FROM shortest_path_DWithin UNION SELECT node_id2 FROM shortest_path_DWithin )  as sub  
			ORDER BY node_id ASC 
		)
		, mapping_lab_to_node AS (
			SELECT art.node_id, rc.cluster_id
			FROM result_clustering AS rc
				LEFT OUTER JOIN artificial_seq AS art ON (rc.seq = art.seq)
		)
		SELECT *
		FROM mapping_lab_to_node as mapl
			LEFT OUTER JOIN patch_node_table AS pn1 USING (node_id )  ;
*/
 
	DROP TABLE IF EXISTS network_clustering ; 
	CREATE TABLE network_clustering AS 
	WITh artificial_seq AS (
		SELECT DISTINCT (row_number() over(order by node_id) -1)::int as seq,   node_id 
			FROM (SELECT node_id1 AS node_id FROM shortest_path_DWithin WHERE tcost >0 UNION SELECT node_id2 FROM shortest_path_DWithin WHERE tcost >0)  as sub  
			ORDER BY node_id ASC 
	)
	, translating_input AS (
		SELECT pn1.seq AS node_id1, pn2.seq AS node_id2, tcost 
		FROM shortest_path_DWithin AS sp
			LEFT OUTER JOIN artificial_seq AS pn1 ON (pn1.node_id = sp.node_id1) 
			LEFT OUTER JOIN artificial_seq AS pn2 ON (pn2.node_id = sp.node_id2) 
	)
	, preparing_input AS (
		SELECT array_agg(node_id1 ORDER BY node_id1,node_id2) as node_id1s
			, array_agg(node_id2 ORDER BY node_id1,node_id2) as node_id2s
			, array_agg(exp(-pow(tcost/38,2) ) ORDER BY node_id1, node_id2) as tcosts --tp get affinity
			-- ,  array_agg( - tcost/285.968150629562  ORDER BY node_id1, node_id2) as tcosts --to get distance
		FROM translating_input
		-- WHERE tcost < 3
	)
	 , result_clustering AS (
	SELECT f.*
	FROM preparing_input,  rc_graph_clustering(node_id1s,node_id2s,tcosts
		,n_clusters := 80
		,eps := 10
		, min_samples := 3
		,preference := 20
		) AS f
	)
	, mapping_lab_to_node AS (
			SELECT art.node_id, rc.cluster_id
			FROM result_clustering AS rc
				LEFT OUTER JOIN artificial_seq AS art ON (rc.seq = art.seq)
	)	
	SELECT *
	FROM mapping_lab_to_node as mapl
		LEFT OUTER JOIN patch_node_table AS pn1 USING (node_id )  ;


-- visualisation and creating a topology	
	DROP TABLE IF EXISTS network_clustering_convexhull ; 
	CREATE TABLE network_clustering_convexhull AS 
	WITH unioned AS (
		SELECT cluster_id, rc_Largest_poly(ST_Buffer(ST_Buffer(ST_Union(ST_Buffer(patch_centroid,8)),-8),2)) as u
			, ST_Collect(patch_centroid) AS centroids
		FROM network_clustering
		GROUP BY cluster_id
	)
	SELECT  cluster_id, ST_Multi(u)::geometry(Multipolygon ,932011 ) AS convex_hull, centroids
	FROM unioned 
	WHERE ST_Area(u) > 50
 

SELECT DropTopology('patch_connectivity_clustering_2') ; 
SELECT CreateTopology('patch_connectivity_clustering_2', 932011, 0.1, hasz:=False); --2 

	SELECT n1.cluster_id AS c1 ,n2.cluster_id AS c2
		, TopoGeo_AddLineString('patch_connectivity_clustering_2'
			, ST_MakeLine(ST_Centroid(n1.centroids), ST_Centroid (n2.centroids)) 
		 ,0.1)
	FROM network_clustering_convexhull AS n1
		,network_clustering_convexhull AS n2
	WHERE ST_DWithin(n1.convex_hull ,n2.convex_hull,1) = TRUE

 SELECT rc_heal_edge_sequentially('patch_connectivity_clustering') ; --note : ùmust be done iteratively ! 





 -------------------------
--  Exporting graph as graphml for external clustering

DROP TABLE IF EXISTS network_clustering ; 
	CREATE TABLE network_clustering AS 

SELECT *
FROM edge_table
LIMIT 10

SELECT *
FROM patch_node_table
LIMIT 10

	WITH artificial_seq AS (
		SELECT DISTINCT (row_number() over(order by node_id) -1)::int as seq,   node_id 
				, round(st_x(centre)::numeric,3) as x 
				, round(st_y(centre)::numeric,3) AS y 
				, round(st_z(centre)::numeric,3) AS z
				FROM patch_node_table , ST_Transform(St_SetSRID(patch_centroid,931008),932011) AS centre
			ORDER BY node_id ASC  
		 --LIMIT 1000
	)
	, node_input AS(
	SELECT array_agg(seq::int ORDER BY seq ASC) AS nodes
		, array_agg(x ORDER BY seq ASC) AS xs
		, array_agg(y ORDER BY seq ASC) AS ys
		, array_agg(z ORDER BY seq ASC) AS zs
		FROM artificial_seq
	)
	 , translating_input AS (
		SELECT pn1.seq AS node_id1, pn2.seq AS node_id2
			,  COALESCE(tcost,0.001)::real AS tcost
		FROM edge_table AS sp
			INNER JOIN artificial_seq AS pn1 ON (pn1.node_id = sp.source) 
			INNER JOIN artificial_seq AS pn2 ON (pn2.node_id = sp.target) 
		 LIMIT 1000
	)
	, preparing_input AS (
		SELECT array_agg(node_id1::int ORDER BY node_id1,node_id2) as node_id1s
			, array_agg(node_id2::int ORDER BY node_id1,node_id2) as node_id2s
			, array_agg(tcost ORDER BY node_id1, node_id2) as tcosts --tp get affinity
			
		FROM translating_input
		-- WHERE tcost < 3
	)
	 
	SELECT f.*
	FROM node_input, preparing_input,   rc_to_graphML ( nodes,xs,ys, zs, node_id1s,node_id2s,tcosts
		 , '/ExportPointCloud/patch_graph.graphml'
		) AS f
	 


 -------------------------