5  Future work & suggestions

5.1 Point of interest (POI)

5.2 Data

Show the code

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")
)

             analysis
1 Vulnerability index
                                                                      problem
1 Results and conclussions based on aggregated data can be limited and biased
                                                                         suggestion
1 Find non-spatial data that can be added using municipality, neighbor or sector ID

5.2.1 Bairros

Show the code

# 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;

5.3 Perforamance

5.3.1 Flood mask

Using the dilate and erode method described in the slide 95/187 (link) could remove small islands that increase the computational costs. A use case of this technique is observed when simplifying coastlines (link)

5.3.2 Buildings

5.3.2.1 Download from Overture

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)

5.3.2.2 Processing in duckdb

Show the code

--- 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');

5.4 Alternative route

Show the code

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 and post

Show the code

---- 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”.

Show the code

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

Show the code

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")  

5.4.1 Porto Alegre

Show the code

--- 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

5.4.2 Canoas

Show the code

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

5.5 Isochrones

5.5.1 Import data

Show the code

# 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
}

5.6 Workflow in R

Show the code

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,])

Show the code

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() 
}
}

Show the code

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))

Show the code

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)

5.7 Visualization

5.8 Tabla A: alternative

Show the code

#### 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)

5.9 Tabla B: centrality

Show the code

---- 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

5.9.1 Data

Show the code

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;

5.9.2 plot

Show the code

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)
  )