future_work.qmd

# Future work & suggestions
## Point of interest (POI)
## Data

```{r}
data.frame(
  analysis = c("Vulnerability index"),
  problem =c("Results and conclussions based on aggregated data can be limited and biased"),
  suggestion =c("Find non-spatial data that can be added using municipality, neighbor or sector ID")
)
```

### Bairros

```{sql}
#| eval: false

# import the data : shp2pgsql -D -I -s 4674 -W "LATIN1" '43SEE250GC_SIR.shp' barrios_brasil | psql -p 25432 -U docker -d gis -h localhost

CREATE TABLE brasil_porto_barrios AS
SELECT 
	nm_bairro,
	st_union(brasil_barrios_4326.geom_4326)
FROM
	brasil_barrios_4326,
	ghs_aoi_porto_alegre
WHERE 
 brasil_barrios_4326.nm_bairro IS NOT NULL
 AND
 	st_intersects(brasil_barrios_4326.geom_4326,
 				ghs_aoi_porto_alegre.geom)
GROUP BY 
	nm_bairro;
```

## Perforamance

### Flood mask

Using the dilate and erode method described in the slide 95/187 ([link](https://cdn2.hubspot.net/hubfs/2283855/PostGIS%20Day%202019%20-%20Overview.pdf)) could remove small islands that increase the computational costs. A use case of this technique is observed when simplifying coastlines ([link](https://github.com/rruiz-s/heigit-gima/wiki/2.-Data-preparation#post-scenario))

### Buildings
#### Download from Overture


```{r}
#| eval: false
#| code-summary: Import buildings and roads from Overture
## Import data
library(overturemapsr)
ghs <- st_read(eisenberg_connection, "ghs_aoi_porto_alegre")
overture_roads <- record_batch_reader(overture_type = 'segment', bbox = sf::st_bbox(ghs))
overture_buildings <- record_batch_reader(overture_type = 'building', bbox = sf::st_bbox(ghs))
## Subset data
overture_roads_subset <- overture_roads |> dplyr::select(c(id, geometry,class,update_time, subtype))
## Export data
### Local file
library(arrow)
library(sfarrow)
sfarrow::st_write_parquet(overture_roads, "overture_roads.parquet")
### Database
library(DBI)
library(nanoarrow)
DBI::dbWriteTable(eisenberg_connection, "overture_roads_subset", overture_roads_subset)
```


#### Processing in duckdb

```{sql}
#| eval: false

--- Using Duckdb
INSTALL spatial;
LOAD spatial;
----- Importing buildings & flood extent from file to duckdb

---- Buildings
CREATE TABLE porto_buildings AS
    SELECT 
        id,
        st_geomfromwkb(geometry) as geom
    FROM 
        'porto_buildings_overture.parquet';
        
--- Floodig      
CREATE TABLE  flood AS
    SELECT
        * 
    FROM 
      st_read('flooding_porto_cleaned.geojson');
---- Filtering building in flood extent
CREATE TABLE flood_building AS 
    SELECT 
      buildings.*
    FROM
      porto_buildings AS buildings,
      flood 
    WHERE
      st_intersects(buildings.geom, flood.geom);
--- Export data
COPY flood_building TO
  'flood_building_v2.gpkg' WITH (FORMAT GDAL, DRIVER 'GPKG');
COPY porto_buildings TO 
   'porto_buildings.gpkg' WITH (FORMAT GDAL, DRIVER 'GPKG');

```


## Alternative route

```{sql}
#| connection: eisenberg_connection
#| eval: false


CREATE TABLE second_path_net AS
WITH first_path AS (
SELECT * FROM pgr_dijkstra(
			'SELECT id,
					 source,
					 target, 
					 cost 
			FROM porto_alegre_net_pre_component_one',
			7249, 7268,
			true))
SELECT 
	* 
FROM 
	porto_alegre_net_pre_component_one
WHERE
	porto_alegre_net_pre_component_one.id NOT IN (select edge from first_path);

SELECT * FROM pgr_dijkstra(
			'SELECT id,
					 source,
					 target, 
					 cost 
			FROM second_path_net',
			7249, 7268,
			true)
-----2


CREATE TABLE pkt_2000 AS
WITH hospital AS (SELECT 
	h.ds_cnes,
	net.id,
	h.geom
FROM hospital_rs_geoportal as h
LEFT JOIN LATERAL
	(SELECT 
		id,
		the_geom
	 FROM
	 porto_alegre_net_pre_component_one_vertices_pgr as net
	 ORDER BY h.geom <-> net.the_geom
	 LIMIT 1) AS net ON TRUE
WHERE 
	h.ds_cnes = 'IRMANDADE DA SANTA CASA DE MISERICORDIA DE PORTO ALEGRE'),
origin AS(SELECT * FROM porto_alegre_net_pre_component_one_vertices_pgr 
WHERE  id = 7268)
SELECT j.path_id, j.node, j.edge, b.cost, b.the_geom FROM pgr_KSP(
						'SELECT id,
					        source,
							target,
							cost 
						FROM porto_alegre_net_pre_component_one',
						(SELECT id FROM hospital),
						(SELECT id FROM origin),
						2000) as j
						left JOIN porto_alegre_net_pre_component_one AS b
				  ON j.edge = b.id;

CREATE TABLE road_without_first AS 
SELECT 
	* 
FROM 
	porto_alegre_net_pre_component_one  
WHERE 
	id NOT IN (SELECT edge FROM pkt_2000 WHERE path_id = 1);
	
	-----3
	WITH paths AS (
  SELECT 
    path_id,
    edge,
    agg_cost,
    SUM(agg_cost) OVER (PARTITION BY path_id ORDER BY seq ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW) AS cum_sum,
    SUM(agg_cost) OVER (PARTITION BY path_id) AS total_cost
  FROM 
    your_table
),
filtered_paths AS (
  SELECT 
    path_id,
    edge,
    agg_cost,
    cum_sum,
    total_cost
  FROM 
    paths
  WHERE 
    cum_sum > 0.25 * total_cost -- remove the first 25%
    AND cum_sum < 0.75 * total_cost -- remove the last 25%
)
SELECT 
  path_id, 
  edge, 
  agg_cost 
FROM 
  filtered_paths;
  
  ----5
  WITH hospital AS (SELECT 
	h.ds_cnes,
	net.id,
	h.geom
FROM hospital_rs_geoportal as h
LEFT JOIN LATERAL
	(SELECT 
		id,
		the_geom
	 FROM
	 porto_alegre_net_pre_component_one_vertices_pgr as net
	 ORDER BY h.geom <-> net.the_geom
	 LIMIT 1) AS net ON TRUE
WHERE 
	h.ds_cnes = 'IRMANDADE DA SANTA CASA DE MISERICORDIA DE PORTO ALEGRE'),
origin AS(SELECT
				* 
		  FROM 
		  		porto_alegre_net_pre_component_one_vertices_pgr 
WHERE  id = 7268),
table_ksp AS (
SELECT j.seq, j.path_id, j.path_seq, j.node, j.edge, j.agg_cost,j.cost,
	b.the_geom, b.target, b.source, b.id FROM pgr_KSP(
						'SELECT id,
					        source,
							target,
							cost 
						FROM porto_alegre_net_pre_component_one',
						(SELECT id FROM hospital),
						(SELECT
							id 
						FROM 
							origin,
							extract_ksp_segments
						WHERE
							origin.id NOT IN (
								SELECT 
									id 
								FROM 
									extract_ksp_segments
								WHERE
									path_id = 1)),
						3) as j
						left JOIN porto_alegre_net_pre_component_one AS b
				  ON j.edge = b.id)
SELECT  * FROM table_ksp;




```

This second one does work without creating two segments of networks. It just increased the costs dinamically. This is inspired by [post](https://gis-ops.com/pgrouting-customize-routing-pleasant-hiking/) and [post](https://gis.stackexchange.com/questions/432741/creating-two-costs-in-pgrouting)

```{sql}
#| eval: false

---- Obtain the shortest path as usual
CREATE TABLE first_path_dijkstra AS
SELECT 
		seq,
		path_seq,
		start_vid,
		end_vid,
		node,
		edge,
		net.the_geom 
FROM  pgr_dijkstra('
				SELECT 
						id,
						source,
						target,
						cost
				FROM porto_alegre_net_largest',
				9372,
				1084) AS path
	LEFT JOIN porto_alegre_net_pre_component_one AS net
	ON path.edge = net.id;
---- detected the shortest path and then increased cost two fold
CREATE TABLE dijkstra_second_route AS
WITH pgr_ksp_table AS (
SELECT 
	*
FROM  pgr_ksp('
				SELECT 
						id,
						source,
						target,
						cost
				FROM porto_alegre_net_largest',
				9372,
				1084,
				3,
				true)), ----node 9372, network = source, target)
second_route AS (
			SELECT 
			net.*,
			CASE 
				WHEN route.node = net.source THEN net.cost * 2
				ELSE net.cost
			END AS cost_updated
			FROM 			
				porto_alegre_net_largest AS net
			LEFT JOIN 
				pgr_ksp_table  AS route 
			ON 	
				route.edge= net.id)	
		SELECT * FROM second_route;

---- Create spatial index	
CREATE INDEX idx_ ON dijkstra_second_route USING gist(the_geom);
CREATE INDEX dijkstra_second_route_source_source ON dijkstra_second_route USING btree(source);
CREATE INDEX dijkstra_second_route_source_target ON dijkstra_second_route USING btree(target);

---- obtain the second shortest path, since the first-path has now higher 
CREATE TABLE second_path_dikstra AS
SELECT 
			seq,
			path_seq,
			start_vid,
			end_vid,
			node,
			edge,
			net.the_geom 
		FROM 
			pgr_Dijkstra('
					SELECT
						id,
						source,
						target,
						cost_updated AS cost
					FROM
						dijkstra_second_route',
						9372,
						1084)  AS path
		LEFT JOIN
				dijkstra_second_route AS net ON
			path.edge = net.id


```

Using this strategy, the following code create 3 alternative routes obtaining the final table "paths_linestring_table".

```{sql}
#| eval: false
CREATE TABLE first_path_6211 AS
SELECT 
		seq,
		path_seq,
		start_vid,
		end_vid,
		node,
		edge,
		net.the_geom,
		1 AS path
FROM  pgr_dijkstra('
				SELECT 
						id,
						source,
						target,
						cost
				FROM porto_alegre_net_largest',
				6211,
				ARRAY(SELECT id FROM hospital_rs_node_v2),
				directed := TRUE) AS path
	LEFT JOIN porto_alegre_net_largest AS net
	ON path.edge = net.id;
---- detected the shortest path and then increased cost two fold

DROP TABLE dijkstra_second_route_6211;
CREATE TABLE dijkstra_second_route_6211 AS
WITH pgr_doubled_dijkstra AS (
SELECT 
	*
FROM  pgr_dijkstra('
				SELECT 
						id,
						source,
						target,
						cost
				FROM porto_alegre_net_largest',
				6211,
				ARRAY(SELECT id FROM hospital_rs_node_v2),
				directed:=true)), ----node 9372, network = source, target)
second_route AS (
			SELECT 
			net.*,
			CASE 
				WHEN route.node = net.source THEN net.cost * 20
				ELSE net.cost
			END AS cost_updated
			FROM 			
				porto_alegre_net_largest AS net
			LEFT JOIN 
				pgr_doubled_dijkstra  AS route 
			ON 	
				route.edge= net.id)	
		SELECT * FROM second_route;

---- Create spatial index	
CREATE INDEX idx_second_6211 ON dijkstra_second_route_6211 USING gist(the_geom);
CREATE INDEX idx_second_6211_source ON dijkstra_second_route_6211 USING btree(source);
CREATE INDEX idx_second_6211_target ON dijkstra_second_route_6211 USING btree(target);

---- obtain the second shortest path, since the first-path has now higher 
DROP TABLE second_path_dikstra_6211;
CREATE TABLE second_path_dikstra_6211 AS
SELECT 
			seq,
			path_seq,
			start_vid,
			end_vid,
			node,
			edge,
			net.the_geom,
			2 AS path
		FROM 
			pgr_Dijkstra('
					SELECT
						id,
						source,
						target,
						cost_updated AS cost
					FROM
						dijkstra_second_route_6211',
					6211,
					ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
		LEFT JOIN
				dijkstra_second_route_6211 AS net ON
			path.edge = net.id;
			
			select * from second_path_dikstra_6211;
		
		
		
----- 3 path
SELECT * FROM dijkstra_second_route_6211;		

DROP TABLE dijkstra_third_route_6211;
CREATE TABLE dijkstra_third_route_6211 AS
WITH pgr_doubled_dijkstra AS (
SELECT 
	*
FROM  pgr_dijkstra('
				SELECT 
						id,
						source,
						target,
						cost_updated AS cost
				FROM dijkstra_second_route_6211',
				6211,
				ARRAY(SELECT id FROM hospital_rs_node_v2),
				directed:=true)), ----node 9372, network = source, target)
second_route AS (
			SELECT 
			net.*,
			CASE 
				WHEN route.node = net.source THEN net.cost_updated * 20
				ELSE net.cost_updated
			END AS cost_updated_nd
			FROM 			
				dijkstra_second_route_6211 AS net
			LEFT JOIN 
				pgr_doubled_dijkstra  AS route 
			ON 	
				route.edge= net.id)	
		SELECT * FROM second_route;
	----

DROP TABLE third_path_dikstra_6211 ;
CREATE TABLE third_path_dikstra_6211 AS
SELECT 
			seq,
			path_seq,
			start_vid,
			end_vid,
			node,
			edge,
			net.the_geom,
			3 AS path
		FROM 
			pgr_Dijkstra('
					SELECT
						id,
						source,
						target,
						cost_updated_nd AS cost
					FROM
						dijkstra_third_route_6211',
					6211,
					ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
		LEFT JOIN
				dijkstra_third_route_6211 AS net ON
			path.edge = net.id;
			
			select * from dijkstra_third_route_6211;	
		

CREATE TABLE three_alternative_paths AS 
SELECT * FROM first_path_6211
UNION
SELECT * FROM second_path_dikstra_6211
UNION
SELECT * FROM third_path_dikstra_6211

CREATE TABLE three_alternative_paths_20 AS 
SELECT * FROM first_path_6211
UNION
SELECT * FROM second_path_dikstra_6211
UNION
SELECT * FROM third_path_dikstra_6211

CREATE TABLE ernesto_three_alternatives AS
SELECT
alternative_paths.*,
hospitals.ds_cnes,
hospitals.cd_cnes
FROM 
three_alternative_paths AS alternative_paths
JOIN hospital_rs_node_v2 AS hospitals
ON hospitals.id = alternative_paths.end_vid
WHERE end_vid = 1743;


CREATE TABLE ernesto_three_alternatives_20 AS
SELECT
alternative_paths.*,
hospitals.ds_cnes,
hospitals.cd_cnes
FROM 
three_alternative_paths_20 AS alternative_paths
JOIN hospital_rs_node_v2 AS hospitals
ON hospitals.id = alternative_paths.end_vid
WHERE end_vid = 1743;


CREATE TABLE paths_linestring_table AS
WITH paths_linestring AS (
SELECT
	start_vid,
	end_vid,
	path,
	st_union(the_geom) AS line_geom
FROM 
	three_alternative_paths_20
GROUP BY path, start_vid, end_vid)
SELECT 
	paths_linestring.*,
	hospitals.ds_cnes,
	round(st_length(line_geom::geography)::numeric/1000,2) AS distance
FROM 
	paths_linestring
JOIN hospital_rs_node_v2 AS hospitals
ON hospitals.id = paths_linestring.end_vid

```


This code created the table to represent these values

```{r}
#| eval: false


alternative_routes <- st_read(eisenberg_connection ,"paths_linestring_table")


library(gtExtras)
alternative_routes |> arrange(distance) |> 
      st_drop_geometry() |>
      select(c(ds_cnes, path, distance)) |> 
      mutate(ds_cnes = as_factor(str_to_title(ds_cnes))) |>
      group_by(ds_cnes) |> 
      mutate(sd = sd(distance)) |> 
    gt(groupname_col="ds_cnes")  

```

### Porto Alegre

```{sql}
#| eval: false

--- Before flooding
CREATE TABLE centroid_cbd_porto AS
SELECT 
	vertices.id,
	vertices.the_geom
FROM 
	central_business_district  AS cbd
LEFT JOIN LATERAL 
	(SELECT
		id,
		the_geom
	 FROM
	 	porto_alegre_net_pre_vertices_pgr AS vertices
	 	ORDER BY 
	 		st_centroid(cbd.st_union) <-> vertices.the_geom
	 		LIMIT 1) AS vertices ON TRUE	
WHERE 
	nm_bairro 
LIKE
	'Centro Histórico';


CREATE TABLE first_path_porto_alegre_cbd AS
SELECT 
        seq,
        path_seq,
        start_vid,
        end_vid,
        node,
        edge,
        net.the_geom,
        1 AS path
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost
                FROM porto_alegre_net_largest',
                (SELECT id FROM centroid_cbd_porto),
                ARRAY(SELECT id FROM hospital_rs_node_v2),
                directed := TRUE) AS path
    LEFT JOIN porto_alegre_net_largest AS net
    ON path.edge = net.id;
    
   
CREATE TABLE dijkstra_second_route_cbd_porto AS
WITH pgr_doubled_dijkstra AS (
SELECT 
    *
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost
                FROM porto_alegre_net_largest',
                 (SELECT id FROM centroid_cbd_porto),
                ARRAY(SELECT id FROM hospital_rs_node_v2),
                directed:=true)), 
second_route AS (
            SELECT 
            net.*,
            CASE 
                WHEN route.node = net.source THEN net.cost * 20
                ELSE net.cost
            END AS cost_updated
            FROM            
                porto_alegre_net_largest AS net
            LEFT JOIN 
                pgr_doubled_dijkstra  AS route 
            ON  
                route.edge= net.id) 
        SELECT * FROM second_route; 
--- Create spatial index   
CREATE INDEX idx_second_cdb_porto ON dijkstra_second_route_cbd_porto USING gist(the_geom);
CREATE INDEX idx_second_cdb_porto_source ON dijkstra_second_route_cbd_porto USING btree(source);
CREATE INDEX idx_second_cdb_porto_target ON dijkstra_second_route_cbd_porto USING btree(target);
  
---- obtain the second shortest path, since the first-path has now higher 
DROP TABLE second_path_dikstra_cdb_porto;
CREATE TABLE second_path_dikstra_cdb_porto AS
SELECT 
            seq,
            path_seq,
            start_vid,
            end_vid,
            node,
            edge,
            net.the_geom,
            2 AS path
        FROM 
            pgr_Dijkstra('
                    SELECT
                        id,
                        source,
                        target,
                        cost_updated AS cost
                    FROM
                        dijkstra_second_route_cbd_porto',
                    (SELECT id FROM centroid_cbd_porto),
                    ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
        LEFT JOIN
                dijkstra_second_route_cbd_porto AS net ON
            path.edge = net.id;
select * from second_path_dikstra_cdb_porto;
           
CREATE TABLE dijkstra_third_route_cdb_porto AS
WITH pgr_doubled_dijkstra AS (
SELECT 
    *
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost_updated AS cost
                FROM dijkstra_second_route_cbd_porto',
                (SELECT id FROM centroid_cbd_porto),
                ARRAY(SELECT id FROM hospital_rs_node_v2),
                directed:=true)), ----node 9372, network = source, target)
second_route AS (
            SELECT 
            net.*,
            CASE 
                WHEN route.node = net.source THEN net.cost_updated * 20
                ELSE net.cost_updated
            END AS cost_updated_nd
            FROM            
                dijkstra_second_route_cbd_porto AS net
            LEFT JOIN 
                pgr_doubled_dijkstra  AS route 
            ON  
                route.edge= net.id) 
        SELECT * FROM second_route;
    ----
       
CREATE TABLE third_path_dikstra_cdb_porto AS
SELECT 
            seq,
            path_seq,
            start_vid,
            end_vid,
            node,
            edge,
            net.the_geom,
            3 AS path
        FROM 
            pgr_Dijkstra('
                    SELECT
                        id,
                        source,
                        target,
                        cost_updated_nd AS cost
                    FROM
                        dijkstra_third_route_cdb_porto',
                    (SELECT id FROM centroid_cbd_porto),
                    ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
        LEFT JOIN
                dijkstra_third_route_cdb_porto AS net ON
            path.edge = net.id;
            
            select * from third_path_dikstra_cdb_porto;    

CREATE TABLE three_alternative_paths_porto AS 
SELECT * FROM first_path_porto_alegre_cbd
UNION
SELECT * FROM second_path_dikstra_cdb_porto
UNION
SELECT * FROM third_path_dikstra_cdb_porto


CREATE TABLE paths_linestring_table_cdb_porto AS
WITH paths_linestring AS (
SELECT
    start_vid,
    end_vid,
    path,
    st_union(the_geom) AS line_geom
FROM 
    three_alternative_paths_porto
GROUP BY path, start_vid, end_vid)
SELECT 
    paths_linestring.*,
    hospitals.ds_cnes,
    round(st_length(line_geom::geography)::numeric/1000,2) AS distance
FROM 
    paths_linestring
JOIN hospitales AS hospitals
ON hospitals.id = paths_linestring.end_vid


SELECT * FROM paths_linestring_table_cdb_porto;

---- after flooding
    
 ---- Postflooding

---- hospitals
 CREATE TABLE hospitals_after_flooding AS 
 SELECT 
 	hospitals.id 
 FROM 
 	hospital_rs_node_v2 AS hospitals,
	flooding_cleaned_porto_union AS flood
WHERE st_disjoint(hospitals.geom_node,
				flood.geom);
	
-----
		
DROP TABLE first_path_porto_alegre_cbd_after;
CREATE TABLE first_path_porto_alegre_cbd_after AS
SELECT 
        seq,
        path_seq,
        start_vid,
        end_vid,
        node,
        edge,
        net.the_geom,
        1 AS path
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost
                FROM  prueba_largest_network_post',
                (SELECT id FROM centroid_cbd_porto),
                ARRAY(SELECT id FROM hospitals_after_flooding),
                directed := TRUE) AS path
    LEFT JOIN prueba_largest_network_post AS net
    ON path.edge = net.id; 
----------------
DROP TABLE dijkstra_second_route_cbd_porto_after;
CREATE TABLE dijkstra_second_route_cbd_porto_after AS
WITH pgr_doubled_dijkstra AS (
SELECT 
    *
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost
                FROM prueba_largest_network_post',
                 (SELECT id FROM centroid_cbd_porto),
                ARRAY(SELECT id FROM hospitals_after_flooding),
                directed:=true)), 
second_route AS (
            SELECT 
            net.*,
            CASE 
                WHEN route.node = net.source THEN net.cost * 20
                ELSE net.cost
            END AS cost_updated
            FROM            
                prueba_largest_network_post AS net
            LEFT JOIN 
                pgr_doubled_dijkstra  AS route 
            ON  
                route.edge= net.id) 
        SELECT * FROM second_route; 
-----
  
SELECT DISTINCT(end_vid) FROM first_path_porto_alegre_cbd_after; ---15 hospitals
SELECT DISTINCT (cd_cnes) FROM hospital_rs_node_v2 ; --- 22 hospital totals

-------------
DROP TABLE second_path_dikstra_cdb_porto_after;	
CREATE TABLE second_path_dikstra_cdb_porto_after AS
SELECT 
            seq,
            path_seq,
            start_vid,
            end_vid,
            node,
            edge,
            net.the_geom,
            2 AS path
        FROM 
            pgr_Dijkstra('
                    SELECT
                        id,
                        source,
                        target,
                        cost_updated AS cost
                    FROM
                        dijkstra_second_route_cbd_porto_after',
                    (SELECT id FROM centroid_cbd_porto),
                    ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
        LEFT JOIN
                dijkstra_second_route_cbd_porto_after AS net ON
            path.edge = net.id;
           

DROP TABLE dijkstra_third_route_cdb_porto_after;
CREATE TABLE dijkstra_third_route_cdb_porto_after AS
WITH pgr_doubled_dijkstra AS (
SELECT 
    *
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost_updated AS cost
                FROM dijkstra_second_route_cbd_porto_after',
                (SELECT id FROM centroid_cbd_porto),
                ARRAY(SELECT id FROM hospital_rs_node_v2),
                directed:=true)), ----node 9372, network = source, target)
second_route AS (
            SELECT 
            net.*,
            CASE 
                WHEN route.node = net.source THEN net.cost_updated * 20
                ELSE net.cost_updated
            END AS cost_updated_nd
            FROM            
                dijkstra_second_route_cbd_porto_after AS net
            LEFT JOIN 
                pgr_doubled_dijkstra  AS route 
            ON  
                route.edge= net.id) 
        SELECT * FROM second_route;
       
-----
DROP TABLE third_path_dikstra_cdb_porto_after;
CREATE TABLE third_path_dikstra_cdb_porto_after AS
SELECT 
            seq,
            path_seq,
            start_vid,
            end_vid,
            node,
            edge,
            net.the_geom,
            3 AS path
        FROM 
            pgr_Dijkstra('
                    SELECT
                        id,
                        source,
                        target,
                        cost_updated_nd AS cost
                    FROM
                        dijkstra_third_route_cdb_porto_after',
                    (SELECT id FROM centroid_cbd_porto),
                    ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
        LEFT JOIN
                dijkstra_third_route_cdb_porto_after AS net ON
            path.edge = net.id;

           
CREATE TABLE three_alternative_paths_porto_after AS 
SELECT * FROM first_path_porto_alegre_cbd_after
UNION
SELECT * FROM second_path_dikstra_cdb_porto_after
UNION
SELECT * FROM third_path_dikstra_cdb_porto_after


CREATE TABLE paths_linestring_table_cdb_porto_after AS
WITH paths_linestring AS (
SELECT
    start_vid,
    end_vid,
    path,
    st_union(the_geom) AS line_geom
FROM 
    three_alternative_paths_porto_after
GROUP BY path, start_vid, end_vid)
SELECT 
    paths_linestring.*,
    hospitals.ds_cnes,
    round(st_length(line_geom::geography)::numeric/1000,2) AS distance
FROM 
    paths_linestring
JOIN hospital_rs_node_v2 AS hospitals
ON hospitals.id = paths_linestring.end_vid

```

### Canoas

```{sql}
#| eval: false



CREATE TABLE centroid_cbd_canoas AS
SELECT 
	vertices.id,
	vertices.the_geom
FROM 
	central_business_district  AS cbd
LEFT JOIN LATERAL 
	(SELECT
		id,
		the_geom
	 FROM
	 	porto_alegre_net_pre_vertices_pgr AS vertices
	 	ORDER BY 
	 		st_centroid(cbd.st_union) <-> vertices.the_geom
	 		LIMIT 1) AS vertices ON TRUE	
WHERE 
	nm_bairro 
LIKE
	'Centro';


CREATE TABLE first_path_canoas_cbd AS
SELECT 
        seq,
        path_seq,
        start_vid,
        end_vid,
        node,
        edge,
        net.the_geom,
        1 AS path
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost
                FROM porto_alegre_net_largest',
                (SELECT id FROM centroid_cbd_canoas),
                ARRAY(SELECT id FROM hospital_rs_node_v2),
                directed := TRUE) AS path
    LEFT JOIN porto_alegre_net_largest AS net
    ON path.edge = net.id;
    

CREATE TABLE dijkstra_second_route_cbd_canoas AS
WITH pgr_doubled_dijkstra AS (
SELECT 
    *
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost
                FROM porto_alegre_net_largest',
                 (SELECT id FROM centroid_cbd_canoas),
                ARRAY(SELECT id FROM hospital_rs_node_v2),
                directed:=true)), 
second_route AS (
            SELECT 
            net.*,
            CASE 
                WHEN route.node = net.source THEN net.cost * 20
                ELSE net.cost
            END AS cost_updated
            FROM            
                porto_alegre_net_largest AS net
            LEFT JOIN 
                pgr_doubled_dijkstra  AS route 
            ON  
                route.edge= net.id) 
        SELECT * FROM second_route; 
--- Create spatial index   
CREATE INDEX idx_second_cdb_porto ON dijkstra_second_route_cbd_porto USING gist(the_geom);
CREATE INDEX idx_second_cdb_porto_source ON dijkstra_second_route_cbd_porto USING btree(source);
CREATE INDEX idx_second_cdb_porto_target ON dijkstra_second_route_cbd_porto USING btree(target);
  
---- obtain the second shortest path, since the first-path has now higher 

CREATE TABLE second_path_dikstra_cdb_canoas AS
SELECT 
            seq,
            path_seq,
            start_vid,
            end_vid,
            node,
            edge,
            net.the_geom,
            2 AS path
        FROM 
            pgr_Dijkstra('
                    SELECT
                        id,
                        source,
                        target,
                        cost_updated AS cost
                    FROM
                        dijkstra_second_route_cbd_canoas',
                    (SELECT id FROM centroid_cbd_canoas),
                    ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
        LEFT JOIN
                dijkstra_second_route_cbd_canoas AS net ON
            path.edge = net.id;
select * from second_path_dikstra_cdb_canoas;
           



CREATE TABLE dijkstra_third_route_cdb_canoas AS
WITH pgr_doubled_dijkstra AS (
SELECT 
    *
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost_updated AS cost
                FROM dijkstra_second_route_cbd_canoas',
                (SELECT id FROM centroid_cbd_canoas),
                ARRAY(SELECT id FROM hospital_rs_node_v2),
                directed:=true)), ----node 9372, network = source, target)
second_route AS (
            SELECT 
            net.*,
            CASE 
                WHEN route.node = net.source THEN net.cost_updated * 20
                ELSE net.cost_updated
            END AS cost_updated_nd
            FROM            
                dijkstra_second_route_cbd_canoas AS net
            LEFT JOIN 
                pgr_doubled_dijkstra  AS route 
            ON  
                route.edge= net.id) 
        SELECT * FROM second_route;
    ----
       
CREATE TABLE third_path_dikstra_cdb_canoas AS
SELECT 
            seq,
            path_seq,
            start_vid,
            end_vid,
            node,
            edge,
            net.the_geom,
            3 AS path
        FROM 
            pgr_Dijkstra('
                    SELECT
                        id,
                        source,
                        target,
                        cost_updated_nd AS cost
                    FROM
                        dijkstra_third_route_cdb_canoas',
                    (SELECT id FROM centroid_cbd_canoas),
                    ARRAY(SELECT id FROM hospital_rs_node_v2))  AS path
        LEFT JOIN
                dijkstra_third_route_cdb_canoas AS net ON
            path.edge = net.id;
            
            select * from third_path_dikstra_cdb_canoas;    

CREATE TABLE three_alternative_paths_canoas AS 
SELECT * FROM first_path_canoas_cbd
UNION
SELECT * FROM second_path_dikstra_cdb_canoas
UNION
SELECT * FROM third_path_dikstra_cdb_canoas


CREATE TABLE paths_linestring_table_cdb_canoas AS
WITH paths_linestring AS (
SELECT
    start_vid,
    end_vid,
    path,
    st_union(the_geom) AS line_geom
FROM 
    three_alternative_paths_canoas
GROUP BY path, start_vid, end_vid)
SELECT 
    paths_linestring.*,
    hospitals.ds_cnes,
    round(st_length(line_geom::geography)::numeric/1000,2) AS distance
FROM 
    paths_linestring
JOIN hospitales AS hospitals
ON hospitals.id = paths_linestring.end_vid



```

## Alternative routes v2

## Isochrones
### Import data

```{r}
#| eval: false

# install.packages("pak")
# pak::pak("GIScience/openrouteservice-r")

library(openrouteservice)

ors_api_key("5b3ce3597851110001cf624887cc526d2c904fb9a2cad4dcd6fe7d31")

poi_hospital <- st_read(eisenberg_connection, "hospital_rs_node_v2")

poi_hospital <- poi_hospital |> arrange(ds_cnes)

### Manually
res_5 <-  ors_isochrones(st_coordinates(poi_hospital$geom_node[1:5]), range = 600, interval = 600, output = "sf") |> mutate(ds_cnes = poi_hospital$ds_cnes[1:5],
                                        cd_cnes = poi_hospital$cd_cnes[1:5],
                                        id = poi_hospital$id[1:5])
res_10 <- ors_isochrones(st_coordinates(poi_hospital$geom_node[6:10]), range = 600, interval = 600, output = "sf") |> mutate(ds_cnes = poi_hospital$ds_cnes[6:10],
                                        cd_cnes = poi_hospital$cd_cnes[6:10],
                                        id = poi_hospital$id[6:10])
res_16 <-  ors_isochrones(st_coordinates(poi_hospital$geom_node[11:15]), range = 600, interval = 600, output = "sf") |>
  mutate(ds_cnes = poi_hospital$ds_cnes[11:15],
        cd_cnes = poi_hospital$cd_cnes[11:15],
                                        id = poi_hospital$id[11:15])
res_22 <-  ors_isochrones(st_coordinates(poi_hospital$geom_node[16:20]), range = 600, interval = 600, output = "sf") |> mutate(ds_cnes = poi_hospital$ds_cnes[16:20],
                                        cd_cnes = poi_hospital$cd_cnes[16:20],
                                        id = poi_hospital$id[16:20])
res_22_24 <-  ors_isochrones(st_coordinates(poi_hospital$geom_node[21:22]), range = 600, interval = 600, output = "sf") |> 
mutate(ds_cnes = poi_hospital$ds_cnes[21:22],
       cd_cnes = poi_hospital$cd_cnes[21:22],
                                        id = poi_hospital$id[21:22])

### 
isochrones <- bind_rows(res_5,res_10, res_16, res_22, res_22_24) |>
  subset(select=c(value, geometry,ds_cnes, cd_cnes,id))

isochrone_v2 <- isochrones |> subset(select=c(group_index, value, geometry)) |> mutate(
  lon=as.numeric(unlist(map(isochrones$center,1))),
  lat=as.numeric(unlist(map(isochrones$center,2))))

DBI::dbWriteTable(eisenberg_connection, "isochrone_600",isochrones)
## 
isochrone_point <- isochrone_v2 |>
                    st_centroid(geometry)
st_join(isochrone_point,
                  poi_hospital, by =st_nearest_features)
##
isochrone_v2 |> mutate(geometry_point = st_as_sf(coords=c("lon","lat")))

poi_hospital_coords <- poi_hospital |> sf::st_drop_geometry() |> 
                    mutate(lon= round(st_coordinates(poi_hospital$geom_node)[,1],3),
lan =st_coordinates(poi_hospital$geom_node)[,2])                            
# sf_use_s2(FALSE)
isochrone_poi <- sf::st_join(sf::st_as_sf(isochrone_v2,poi_hospital_coords), join=st_intersects)

left_join(isochrone_v2, poi_hospital_coords, by="lon")


DBI::dbWriteTable(eisenberg_connection,  "isochrone_v2",isochrone_v2)

### Loop - it does not work
for (i in 1:length(poi_hospital$geom_node)){
  print(i)
   res <- ors_isochrones(
     st_coordinates(poi_hospital$geom_node[i,]),
                     range=1800,
                     interval=600,
                     output="sf")
   isochrone_results[[i]] <- res
}

```

## Workflow in R

```{r}
#| eval: false

library(sfnetworks)
net <- sf::read_sf(eisenberg_connection,
            "porto_alegre_net_largest") |>
  sfnetworks::as_sfnetwork() 

from <- sf::read_sf(eisenberg_connection,
                    "points_cardiologia") |> st_as_sf()
to <- sf::read_sf(eisenberg_connection,
                    "hospital_rs_node_v2") |> st_as_sf()

first_destinations <- st_network_paths(net, from =from[1,], to = to[20,])

```

```{r}
#| eval: false

library(DBI)
library(glue)
library(RPostgres)

eisenberg_connection <- DBI::dbConnect(RPostgres::Postgres(),
                          user= "docker",
                          password = "docker",
                          host = "localhost",
                          dbname="gis",
                          port = 25432)

## Test 2
to <- sf::read_sf(eisenberg_connection,
                  "hospital_rs_node_v2") |> st_as_sf()


table_end_v2 <- glue_sql("CREATE TABLE table_end_v2
                (start_vid integer,
                end_vid integer,
                path smallint,
                length float,
                cost float,
                the_geom geometry)")
table_end_db <- dbGetQuery(conn=eisenberg_connection,table_end_v2 )

isochrones_sampling_points <-st_read(eisenberg_connection, "isochrones_snapped")
df_isochrone_points  <- isochrones_sampling_points 

DBI::dbWriteTable(eisenberg_connection,"isochrones_sampling_points",isochrones_sampling_points)
# we have the id of the hospital, then we nned all the points that are in the isochrone of the hospital. How do we get this?


idx_hospital <- 15
for (idx_hospital in 1:length(unique(isochrones_sampling_points$cd_cnes))){
  tic(glue("hospital run {idx_hospital}/{length(unique(isochrones_sampling_points$cd_cnes))}"))
      code_hospital<-unique(isochrones_sampling_points$cd_cnes)[idx_hospital]
      points_for_isochrone <- filter(isochrones_sampling_points,
                                    cd_cnes == code_hospital)
      alternative_routes(
        df_isochrone_points=points_for_isochrone)
  toc()
  }

#| eval: false

#| eval: false


#idx <- 1



tic.clear()


#| eval: false


library(tictoc)
## data= isochrones_sampling_points
alternative_routes <- function(df_isochrone_points){
for (idx in 1:length(df_isochrone_points$id)) {
  tic(glue("run {idx}/{length(df_isochrone_points$id)}")) 
  hospital_name <- df_isochrone_points$cd_cnes[idx]
  origin_id <- df_isochrone_points$id[idx]
  
  table_name_1 <- glue("first_path_point_{as.character(idx)}")
  table_name_2 <- glue("second_path_point_{as.character(idx)}")
  table_name_3 <- glue("second_route_point_{as.character(idx)}")
  table_name_4 <- glue("third_path_point_{as.character(idx)}")
  table_name_5 <- glue("third_route_point_{as.character(idx)}")
  table_name_6 <- glue("three_alternative_routes_point_{as.character(idx)}")
  table_name_7 <- glue("table_end_point_{as.character(idx)}")
  
  
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`table_name_1`}", .con=eisenberg_connection),
            conn=eisenberg_connection)
  
  # Assuming you have an open PostgreSQL connection in con
  query_1 <- glue_sql("
CREATE TABLE {`table_name_1`} AS
SELECT 
        seq,
        path_seq,
        start_vid,
        end_vid,
        net.the_geom,
        node,
        edge,
        1 AS path
  FROM pgr_dijkstra('
        SELECT
            id,
            source,
            target,
            cost
        FROM porto_alegre_net_largest',
                    ARRAY(SELECT id FROM isochrones_snapped WHERE id ={origin_id}),
                    ARRAY(SELECT id AS end_id FROM hospital_rs_node_v2 WHERE cd_cnes = {hospital_name}),
                    directed:=TRUE) AS path
LEFT JOIN porto_alegre_net_largest AS net
ON path.edge = net.id", .con = eisenberg_connection)
  # Send the query
  result_q1 <- dbGetQuery(conn= eisenberg_connection, query_1)
  # dijkstra_second_route
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`table_name_2`}", .con=eisenberg_connection),
            conn=eisenberg_connection)
  query_2 <- glue_sql("
CREATE TABLE {`table_name_2`} AS
WITH pgr_doubled_dijkstra AS (
    SELECT 
        *
    FROM  pgr_dijkstra('
            SELECT 
                id,
                source,
                target,
                cost
            FROM porto_alegre_net_largest',
            ARRAY(SELECT id FROM isochrones_snapped WHERE id = {origin_id}),
            ARRAY(SELECT id FROM hospital_rs_node_v2 WHERE cd_cnes = {hospital_name}),
            directed := TRUE  -- This clause should be inside the pgr_dijkstra function
        )
),  -- Corrected closing of pgr_dijkstra CTE
second_route AS (
    SELECT 
        net.*,
        CASE 
            WHEN route.node = net.source THEN net.cost * 20
            ELSE net.cost
        END AS cost_updated
    FROM            
        porto_alegre_net_largest AS net
    LEFT JOIN 
        pgr_doubled_dijkstra  AS route 
    ON  
        route.edge = net.id
) 
SELECT * 
FROM second_route", .con = eisenberg_connection)
  
  result_2 <- dbGetQuery(conn= eisenberg_connection, query_2)
  
  ## Calculate teh second route
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`table_name_3`}", .con=eisenberg_connection),
            conn=eisenberg_connection)
  
  query_3 <- glue_sql("
CREATE TABLE {`table_name_3`} AS
SELECT 
    seq,
    path_seq,
    start_vid,
    end_vid,
    node,
    edge,
    net.the_geom,
    2 AS path
FROM 
    pgr_dijkstra('
        SELECT
            id,
            source,
            target,
            cost_updated AS cost
        FROM
            {`table_name_2`}',
        ARRAY(SELECT id FROM isochrones_snapped WHERE id = {origin_id}),
        ARRAY(SELECT id FROM hospital_rs_node_v2 WHERE cd_cnes = {hospital_name}),
        directed:=TRUE
    ) AS path
LEFT JOIN
    {`table_name_2`} AS net 
    ON path.edge = net.id;
", .con = eisenberg_connection)
  
  # Run the CREATE TABLE query
  dbExecute(conn = eisenberg_connection, query_3)
  select_query <- glue_sql("SELECT * FROM {`table_name_3`} ;",.con=eisenberg_connection) 
  
  # Fetch the data
  result_3 <- dbGetQuery(conn = eisenberg_connection, select_query)
  
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`table_name_4`}", .con=eisenberg_connection),
            conn=eisenberg_connection)
  # Third route
  query_4 <- glue_sql("
CREATE TABLE {`table_name_4`} AS
WITH pgr_doubled_dijkstra AS (
SELECT 
    *
FROM  pgr_dijkstra('
                SELECT 
                        id,
                        source,
                        target,
                        cost_updated AS cost
                FROM {`table_name_2`}',
                ARRAY(SELECT id FROM isochrones_snapped WHERE id = {origin_id}),
        ARRAY(SELECT id FROM hospital_rs_node_v2 WHERE cd_cnes = {hospital_name}),
                directed:=TRUE)), ----node 9372, network = source, target)
second_route AS (
            SELECT 
            net.*,
            CASE 
                WHEN route.node = net.source THEN net.cost_updated * 20
                ELSE net.cost_updated
            END AS cost_updated_nd
            FROM            
                {`table_name_2`} AS net
            LEFT JOIN 
                pgr_doubled_dijkstra  AS route 
            ON  
                route.edge= net.id) 
        SELECT * FROM second_route", .con=eisenberg_connection);
  
  result_4 <- dbGetQuery(conn = eisenberg_connection, query_4)
  
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`table_name_5`}", .con=eisenberg_connection),
            conn=eisenberg_connection)
  # third route
  query_5 <- glue_sql("
CREATE TABLE {`table_name_5`} AS
SELECT 
            seq,
            path_seq,
            start_vid,
            end_vid,
            node,
            edge,
            net.the_geom,
            3 AS path
        FROM 
            pgr_Dijkstra('
                    SELECT
                        id,
                        source,
                        target,
                        cost_updated_nd AS cost
                    FROM
                        {`table_name_4`}',
                   ARRAY(SELECT id FROM isochrones_snapped WHERE id = {origin_id}),
        ARRAY(SELECT id FROM hospital_rs_node_v2 WHERE cd_cnes = {hospital_name}),
        directed:=TRUE)  AS path
        LEFT JOIN
                {`table_name_4`} AS net ON
            path.edge = net.id;", .con= eisenberg_connection)
  result <- dbGetQuery(conn = eisenberg_connection, query_5)
  
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`table_name_6`}", .con=eisenberg_connection),
            conn=eisenberg_connection)  
  alterantive_point_1 <- glue_sql("
CREATE TABLE {`table_name_6`} AS 
SELECT start_vid, end_vid,path, the_geom, edge FROM {`table_name_1`}
UNION
SELECT start_vid, end_vid,path, the_geom, edge FROM {`table_name_3`}
UNION
SELECT start_vid, end_vid,path, the_geom, edge FROM {`table_name_5`}", .con = eisenberg_connection)
  
  result_6 <- dbGetQuery(conn = eisenberg_connection, alterantive_point_1)
  
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`table_name_7`}", .con=eisenberg_connection),
            conn=eisenberg_connection)  
  table_end_point_1 <- glue_sql("
CREATE TABLE {`table_name_7`} AS
SELECT
            alternative.*,
            ceil(st_length(alternative.the_geom::geography)/1000) As length,
            original.cost AS cost
          FROM
            {`table_name_6`} AS alternative
          LEFT JOIN
           porto_alegre_net_largest  AS original
          ON
          alternative.edge = original.id
          WHERE
            edge!= -1", .con=eisenberg_connection)
  result_7 <- dbGetQuery(conn=eisenberg_connection,table_end_point_1 )
  
  
  calculations <- glue_sql("
  INSERT INTO table_end_v2 (start_vid, end_vid, path, length, cost, the_geom)
        SELECT 
          start_vid,
          end_vid,
          path,
          sum(length) AS length,
          sum(cost) AS cost,
          st_union(the_geom) as the_geom
        FROM 
            {`table_name_7`} 
        GROUP BY
        start_vid,
        end_vid,
        path
        ORDER BY path ASC", .con=eisenberg_connection)
  
  dbExecute(conn=eisenberg_connection, calculations)
  toc() 
}
}

```


```{r}
#| eval: false

result |>
  dplyr::select(c(start_vid, end_vid,path, length, cost)) |>
  group_by(start_vid,end_vid,path) |>
  summarise(length=sum(length),
            cost=sum(cost))
```


```{sql}
#| connection: eisenberg_connection
#| eval: false

SELECT
                *
        FROM pgr_dijkstra('
        					SELECT
                        id,
                        source,
                        target,
                        cost
                  FROM porto_alegre_net_largest',
               ARRAY(SELECT id FROM points_cardiologia LIMIT 1),
               ARRAY(SELECT id AS end_id FROM hospital_rs_node_v2 WHERE ds_cnes = 'HOSPITAL SAO CAMILOESTEIO'),
               directed:=TRUE)


```



## Visualization
## Tabla A: alternative

```{r}
#| eval: false

#### alternative routes before
alternative_routes <- st_read(eisenberg_connection,                         "paths_linestring_table")
alternative_routes$momment <- "pre"
#### alternative routes after
alternative_routes_after <- st_read(eisenberg_connection,       "paths_linestring_table_cdb_porto_after")
alternative_routes_after$momment <- "post"

##### binding the two
alternative_routes_both <- alternative_routes |> 
  bind_rows(alternative_routes_after)
#### Creating table for before flooding
table_before |> arrange(ds_cnes)  |>  write_csv("table_before_sd.csv")
table_before <- alternative_routes |> 
  rename(distance_before = distance) |> 
  arrange(distance_before) |> 
  st_drop_geometry() |>
  select(c(ds_cnes, path, distance_before)) |> 
  mutate(ds_cnes = as_factor(str_to_title(ds_cnes))) |>
  group_by(ds_cnes) |> 
  mutate(sd_distance_before = sd(distance_before))
##### Creating table for after flooding
table_after <- alternative_routes_after |> arrange(distance) |> 
  st_drop_geometry() |>
  rename(distance_after = distance) |> 
  select(c(ds_cnes, path, distance_after)) |> 
  mutate(ds_cnes = as_factor(str_to_title(ds_cnes))) |>
  group_by(ds_cnes) |> 
  mutate(sd_distance_after = sd(distance_after))
table_after |> arrange(ds_cnes)  |> write_csv("table_after_sd.csv")

#### Creating sd table
tabla_alternative_roads_sd <- alternative_routes_both |>
  sf::st_drop_geometry() |> 
  mutate(ds_cnes = as_factor(str_to_title(ds_cnes)))  |> 
  pivot_wider(names_from = momment,
              values_from=distance,
              names_glue = "{momment}_distance") |> 
  group_by(end_vid,ds_cnes) |>
  summarise(sd_distance_pre = sd(pre_distance, na.rm=TRUE),
            sd_distance_post = sd(post_distance, na.rm =TRUE)) |>
  mutate(change_sd = (sd_distance_post-sd_distance_pre)/sd_distance_pre)

```

## Tabla B: centrality

```{sql}
#| eval: false

---- bairros

DROP TABLE bairro_porto_cleaned ;
create table bairro_porto_cleaned AS
select distinct on (cd_geocodb)
		cd_geocodb,
		tipo,
		nm_bairro,
		st_union(st_transform(geom,4326)) AS the_geom
from 
	barrios_brasil 
where nm_bairro IS NOT NULL
GROUP BY
	cd_geocodb, tipo, nm_bairro;


---- subdistricts cd_geocods 
create table subdistricts_porto_cleaned AS
select distinct on (cd_geocods)
		cd_geocods,
		tipo,
		nm_bairro,
		nm_subdist,
		st_union(st_transform(geom,4326)) AS the_geom
from 
	barrios_brasil 
GROUP BY
	cd_geocods, tipo, nm_bairro,nm_subdist;
---- districts
create table districts_porto_cleaned AS
select distinct on (cd_geocodd)
		cd_geocodd,
		tipo,
		nm_distrit,
		st_union(st_transform(geom,4326)) AS the_geom
from 
	barrios_brasil 
GROUP BY
	cd_geocodd, tipo, nm_distrit;

---

SELECT * FROM bairro_porto_cleaned;
SELECT * FROM barrios_brasil;
SELECT * FROM districts_porto_cleaned;
----- Municipalities
create table municipalities_porto_cleaned AS
select distinct on (cd_geocodm)
		cd_geocodm,
		tipo,
		nm_municip,
		st_union(st_transform(geom,4326)) AS the_geom
from 
	barrios_brasil 
GROUP BY
	cd_geocodm, tipo, nm_municip;
----- centrality para barrios
CREATE TABLE table_barrios_mean AS
WITH centrality_barrios AS (
SELECT 
	cd_geocodb,
	bairros.nm_bairro,
	bairros.tipo,
	centrality.*,
	bairros.the_geom AS polygon
FROM 
	centrality_weighted_100_bidirect_cleaned AS centrality	
LEFT JOIN bairro_porto_cleaned AS bairros
ON st_intersects(centrality.the_geom,
					bairros.the_geom))
SELECT 
	cd_geocodb AS bairro_id,
	nm_bairro AS bairro_name,
	tipo AS bairro_type,
	sum(st_length(the_geom)) AS total_length,
	avg(centrality) AS mean_centrality,
	percentile_disc(0.5) WITHIN GROUP (ORDER BY centrality) AS median,
	polygon
FROM centrality_barrios	
WHERE nm_bairro IS NOT NULL
GROUP BY cd_geocodb, nm_bairro, tipo, polygon


```


### Data 

```{sql}
#| connection: eisenberg_connection
#| eval: false

CREATE TABLE centrality_both_sankey  AS
WITH centrality_both_sankey AS (
SELECT  centrality_pre.bidirectid,
        centrality_pre.centrality AS centrality_pre,
        centrality_post.centrality AS centrality_post,
        centrality_pre.the_geom AS pre_geom,
        centrality_post.the_geom AS post_geom
FROM centrality_weighted_100_bidirect_cleaned AS centrality_pre
FULL OUTER JOIN centrality_weighted_100_bidirect_cleaned_post AS centrality_post
  ON centrality_pre.bidirectid = centrality_post.bidirectid)
SELECT 
    bidirectid,
    centrality_pre,
    centrality_post,
    COALESCE(ST_Length(pre_geom::geography)/100, ST_Length(post_geom::geography)/1000) AS geom_length
FROM 
    centrality_both_sankey;

```

### plot

```{r}
#| eval: false
#| message: false
#| warning: false

library(tidyverse)
library(ggsankey)

### Flow Chart
centrality_both_sankey <- DBI::dbReadTable(eisenberg_connection, "centrality_both_sankey")

df_centrality_both_sankey <- centrality_both_sankey |> 
  mutate(pre_scenario=
           as.factor(case_when(
             centrality_pre == 0 ~ "[0]",
             centrality_pre > 0 & centrality_pre <= 146 ~ "[1-146]",
             centrality_pre >= 146 & centrality_pre <=345 ~ "[146-345]",
             centrality_pre >= 345 & centrality_pre <= 633  ~ "[345-633]",
             centrality_pre >= 633 ~ "[633-1644]",
             TRUE ~ "[0]")),
         post_scenario=
           as.factor(case_when(
             centrality_post == 0 ~ "[0]",
             centrality_post > 0 & centrality_post <= 146 ~ "[1-146]",
             centrality_post >= 146 & centrality_post <=345 ~ "[146-345]",
             centrality_post >= 345 & centrality_post <= 633  ~ "[345-633]",
             centrality_post >= 633 ~ "[633-1644]",
             TRUE ~ "[0]"))) 

# HeiGIt --- Create a long format dataset for ggsankey
df_centrality_both_sankey_long <- df_centrality_both_sankey |> 
  mutate(
    pre_scenario = as.character(pre_scenario),
    post_scenario = as.character(post_scenario)) |> 
  ggsankey::make_long(5:6) 

# Create an initial Sankey diagram (not used in final plot, but used to extract data)
mycolour <- c("#edf8fb","#b3cde3","#8c96c6","#8856a7","#810f7c")
names(mycolour) <- c(unique(df_centrality_both_sankey_long$node)[3],
                     unique(df_centrality_both_sankey_long$node)[1],
                     unique(df_centrality_both_sankey_long$node)[2],
                     unique(df_centrality_both_sankey_long$node)[4],
                     unique(df_centrality_both_sankey_long$node)[5])
auxiliar_data <-df_centrality_both_sankey_long |> group_by(node, x) |> count()
df_centrality_both_sankey_long$node <- factor(df_centrality_both_sankey_long$node, levels=c("[0]",
                                                                                            "[1-146]",
                                                                                            "[146-345]",
                                                                                            "[345-633]",
                                                                                            "[633-1644]")) 
p <- ggplot(df_centrality_both_sankey_long, aes(x = x,
                                                next_x = next_x,
                                                node = node,
                                                next_node = next_node,
                                                fill = factor(node),
                                                label = node)) +
  geom_sankey(flow.alpha = 0.6) +
  geom_sankey_text() +
  scale_fill_manual(values = mycolour) + 
  geom_sankey_label( fill="white", alpha=.7) +
  theme_minimal() +
  theme(
    legend.position = "none",
    plot.background = element_rect(fill = "grey99", color = NA)
  )

# Extract flow, node, and label data from the Sankey diagram
moves_flows <- layer_data(p, 1) %>% 
  mutate(height = flow_end_ymax - flow_end_ymin)

moves_nodes <- layer_data(p, 2)

moves_labels <- layer_data(p, 3)

f1 <- "Graphik"
f1b <- "Graphik Compact"
f2 <- "Publico Headline"

# Apply factor levels to ensure matching with mycolour
moves_nodes$label <- factor(moves_nodes$label, levels = names(mycolour))
moves_flows$label <- factor(moves_flows$label, levels = names(mycolour))

p1 <- ggplot() +
  # Add flow polygons, highlighting the most common moves and using the correct fill mapping
  geom_polygon(data = moves_flows, aes(x, y, group = group, fill = label)) +
  # Add node rectangles with color based on the label
  geom_rect(data = moves_nodes, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = "white"), fill="white") +
  # Manually scale the fill using the mycolour palette
  scale_fill_manual(values = mycolour) +
  # Add labels for the most common moves
  geom_text(data = moves_labels,
            aes(x, y, label = label),
            family = f1b,
            color = "black",
            size = 6) +
  labs(
    title = "Changes on Connectivity",
    subtitle = "Porto Alegre Settlement",
    caption = "Source: OpenStreetMap · Graphic: Ricardo Ruiz Sánchez inspired by Georgios Karamanis"
  ) +
  theme_void(base_family = f1) +
  theme(
    legend.position = "none",
    plot.background = element_rect(fill = "white", color = NA),
    plot.margin = margin(10, 10, 10, 10),
    plot.title = element_text(size = 30, face = "bold", family = f2),
    plot.subtitle = element_text(size = 15, lineheight = 1, color = "#676162"),
    plot.caption = element_text(hjust = 0, size = 12, "#676162")
  )

ggsave("plot_flow_chart.jpg", plot= p1, dpi= 300, width= 10, height= 15)

## Length roads
### pre-length
df_length_pre <- centrality_both_sankey |> select(bidirectid, centrality_pre, geom_length) |>
  rename(centrality = centrality_pre) |> 
  mutate(
  category_centrality=
    as.factor(case_when(
      centrality == 0 ~ "[0]",
      centrality > 0 & centrality <= 146 ~ "[1-146]",
      centrality >= 146 & centrality <=345 ~ "[146-345]",
      centrality >= 345 & centrality <= 633  ~ "[345-633]",
      centrality >= 633 ~ "[633-1644]",
      TRUE ~ "[0]"))) |> 
  dplyr::select(c(geom_length, category_centrality)) |> 
  group_by(category_centrality) |>
  summarise(total_lenth=sum(geom_length))
df_length_pre$event <- "Road length (%)"
df_length_pre$pct <- as.numeric(round(df_length_pre$total_lenth/(sum(df_length_pre$total_lenth)),4))

## pre-length-plot
ggplot(df_length_pre, aes(x = event, y = pct, fill = category_centrality)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = mycolour) +
  #theme_minimal() +
  ggthemes::theme_clean() +
  labs(x="", y="") +
  theme(
    legend.position = "none",
    legend.text = element_text(size = 7),
    legend.title = element_text(size = 7),
    legend.background = element_blank(),
    plot.background = element_blank(),
    plot.caption = element_text(face = "italic",
                                size = 8,
                                color = "grey30"
    )
  ) +
  guides (
    fill = guide_legend(nrow =2, byrow = TRUE)
  ) 

## Centralities values
df_barplots_sankey <- df_centrality_both_sankey |>
  pivot_longer(cols=c("post_scenario","pre_scenario"),
               names_to ="event",
               values_to = "quartile")
centrality_both_longer <- centrality_both |> pivot_longer(cols=c(pre_scenario, post_scenario),
                                names_to = "quartile",
                                values_to = "category") 
df_barplots_sankey <- centrality_both_longer |> 
  dplyr::select(c(centrality, quartile, category)) |>
  group_by(quartile, category) |> summarise(centrality=sum(centrality)) |>
  mutate(total_perc = centrality/sum(centrality))

ggplot(df_barplots_sankey, aes(x = quartile, y = total_perc, fill = category)) +
  geom_col() +
  scale_x_discrete(labels=c("centrality_post" = "Post-flooding", "centrality_pre" = "Pre-flooding")) +
  coord_flip() +
  scale_fill_manual(values = mycolour) +
  labs(
    x = "Event",
    y = "Centrality (%)",
    title = "Share of centrality per range in the two scenarios"
  ) +
  #theme_minimal() +
  ggthemes::theme_clean() +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 7),
    legend.title = element_text(size = 7),
    legend.background = element_blank(),
    plot.background = element_blank(),
    plot.caption = element_text(face = "italic",
                                size = 8,
                                color = "grey30"
    )
  ) +
  guides (
    fill = guide_legend(nrow = 1, byrow = TRUE)
  ) 

```