6 Calculate connectivity

6.1 Research question 1

How is road connectivity affected in pre-disaster and post-disaster scenarios based on edge betweenness centrality?

6.1.1 Core Metropolitan scale

Figure 6.1 shows an interactive view of the connectivity before and after the flooding event in the Core Metropolitan area. On the left side, the urban arteries such as BR-290 freeway¹ is displayed in purple previous to the flood event. On the right side, the flood reduced the connectivity of the urban arteries. Municipalities such as Nova Santa Rita were isolated from the core metropolitan area of Porto Alegre.

Figure 6.1: A comparison of the centrality before and after the floods in the core metropolitan (CM) area

How to visualise CM connectivity

library(mapview)
library(leaflet)
library(leafpop)

# color: https://gis.stackexchange.com/questions/312027/r-and-mapview-how-can-i-chose-different-colors-for-different-layers

### For pre-event: 1644, 468, 142
pal <- c("#b3cde2","#8a95c5","#8854a7","#7a0f76")
centrality_pre$centrality_fct <- cut(centrality_pre$centrality,
                                     breaks=c(0,725,2114,4870,11205),
                                     labels =c("low [1-726]","intermediate [726-2114]","high [2114-4870]","urban arteries [4780-1125]"),
                                     include.lowest= TRUE,
                                     right =FALSE)
#### natural breaks:  81, 230, 582
centrality_post$centrality_fct <- cut(centrality_post$centrality,
                                      breaks=c(0,725,2114,4870,11205),
                                      labels =c("low [1-726]","intermediate [726-2114]","high [2114-4870]","urban arteries [4780-1125]"),
                                      include.lowest= TRUE,
                                      right =FALSE)
### 
library(scales)
centrality_pre$scaled_lwd <- rescale(centrality_pre$centrality, to = c(0.5, 8)) 
centrality_pre_map <- mapview::mapview(centrality_pre,
                                       zcol="centrality_fct",
                                       lwd = centrality_pre$scaled_lwd,
                                       col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                       col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                       map.type=c("OpenStreetMap","CartoDB.DarkMatter"),
                                       layer.name ="Centrality Pre-Event",
                                       popup=leafpop::popupTable(centrality_pre, 
                                                        zcol=c("centrality","bidirectId"))) 
centrality_post_map <- mapview::mapview(centrality_post,
                                        zcol="centrality_fct",
                                        col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                        col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                        map.type=c("OpenStreetMap","CartoDB.DarkMatter"),
                                        lwd = centrality_pre$scaled_lwd,
                                        layer.name ="Centrality Post-Event",
                                        popup=leafpop::popupTable(centrality_post, 
                                                         zcol=c("centrality","bidirectId"))) 

centrality_pre_map | centrality_post_map + mapview::mapview(flood_extent,
                                                   col.regions="#dbc3b0",
                                                   map.type=c("OpenStreetMap","CartoDB.DarkMatter"),
                                                   alpha.regions= 0.5,
                                                   layer.name="Flooding layer")

We used pgrouting to calculate the centrality by counting the least-cost paths using pgr_dijkstra between origin and destinations (OD). The code [ECR 10] is based on Daniel i. Patterson post ² and Matt (Forrest, Mitchell, and Mitchell 2023) The line 7 executed the function pgr_dijkstra from pgrouting creating a table of the shortest paths from the 300 OD contained in the table od_300_snapped_v2. Later, the line 18 added the spatial geometry from the road network saved in the table ghs_aoi_net_largest. A higher number of OD returned the error in handling memory allocation ³, so we limited to 300 origins and destinations. However, the visual inspection of the spatial connectivity patterns were very similar between 100 and 300 OD revealing the BR-290 freeway as the main urban arteria in both cases.

[ACR 10] How to calculate edge betweenness centrality in the Core Metropolitan Area

EXPLAIN ANALYZE
CREATE TABLE centrality_ghs_aoi_pre AS
SELECT
    b.id,
    b.geom,
    count(geom) as centrality
FROM pgr_dijkstra(
    'SELECT 
        id AS id,
        "fromId" AS source,
        "toId" AS target,
        weight as cost
    FROM 
      ghs_aoi_net_largest',
    ARRAY(SELECT snapped_id AS start_id FROM od_300_snapped_v2),
    ARRAY(SELECT snapped_id AS end_id FROM od_300_snapped_v2),
        directed:=TRUE) j
    LEFT JOIN 
      ghs_aoi_net_largest AS b
    ON
      j.edge=b.id
    GROUP BY
        b.id, b.geom;        
        
--- Time: 80542.531 ms

The code chunks from ACR11 to ACR13 handles a characteristic from exporting the ORS road network into R, the bidirectid variable.

[ACR 11] How to add a bidirectId when not present

CREATE TABLE centrality_ghs_aoi_bidirect_pre AS
SELECT 
    t1.*,
    t2."bdrctId"
FROM
    centrality_ghs_aoi_pre AS t1
JOIN
    ghs_aoi_net_largest AS t2
ON
t1.id = t2.id;
--- Time: 164.685 ms

select max("bdrctId") from centrality_ghs_aoi_bidirect_pre ; ---99527
create sequence bidirectid_centrality_ghs_aoi start 99527;
update centrality_ghs_aoi_bidirect_pre
set "bdrctId" = nextval('bidirectid_centrality_ghs_aoi')
where "bdrctId" is null ;
--- checking that bdrctId is not null
select count(*)
from centrality_ghs_aoi_bidirect_pre
where "bdrctId" is null; ---0

We summed up the centrality values of the roads with same bidirectid.

[ACR 12] How to use bidirectid to calculcalte centrality

create table centrality_ghs_aoi_bidirect_pre_group as
select
"bdrctId",
sum(centrality) as centrality
from centrality_ghs_aoi_bidirect_pre
group by "bdrctId";
--- Time: 65.993 ms
create table centrality_ghs_aoi_bidirect_pre_group_id as
select t1.*, t2.id
from centrality_ghs_aoi_bidirect_pre_group t1
join centrality_ghs_aoi_bidirect_pre t2
on t1."bdrctId" = t2."bdrctId";

[ACR 13] How to calculate final edge betweenness centrality in the Core Metropolitan Area

create table centrality_ghs_aoi_bidirect_pre_cleaned as
select
    t1.*,
    t2.geom,
    t2."toId",
    t2."fromId",
    t2.weight,
    t2."undrctI"
from
    centrality_ghs_aoi_bidirect_pre_group_id t1
join
    ghs_aoi_net_largest t2
on
t1.id = t2.id;
--- Time: 565.951 ms

This process is repeated again changing the table containing the road network ghs_aoi_net_largest to ghs_aoi_net_largest_post. While the table ghs_aoi_net_largest represents the road network of the AoI before the flood, the table ghs_aoi_net_largest_post shows the differences on the road network after the flood.

6.1.2 Intracity scale

Figure 6.2 shows the changes on the centrality when we limited the AoI to each municipality, in this case the municipality of Canoas. The flood hit the west side of Canoas affecting the Mathias Velho dam ⁴ causing the loss of urban arterias such as Rua Florianópolis.

Figure 6.2: A comparison of the centrality before and after the floods in the core metropolitan area

How to visualiza a comparison of the centrality before and after the floods in the intracity (IC) area

library(leaflet)
library(mapview)
library(leafpop)

### For pre-event: 1644, 468, 142
pal <- c("#b3cde2","#8a95c5","#8854a7","#7a0f76")
local_connectivity_canoas <- filter(local_connectivity_pre, shapename=="Canoas")
local_connectivity_canoas_post <- filter(local_connectivity_post, shapename=="Canoas")
library(BAMMtools)
natural_jenks_breaks_cmc_pre_canoas <- BAMMtools::getJenksBreaks(local_connectivity_canoas$centrality, k=5)

local_connectivity_canoas$centrality_fct <- cut(local_connectivity_canoas$centrality,
                                     breaks=c(
                                              natural_jenks_breaks_cmc_pre_canoas[1],
                                              natural_jenks_breaks_cmc_pre_canoas[2],
                                              natural_jenks_breaks_cmc_pre_canoas[3],
                                              natural_jenks_breaks_cmc_pre_canoas[4],
                                              natural_jenks_breaks_cmc_pre_canoas[5]),
                                     labels =c("low [1-533]","intermediate [533-1510]","high [1510-2899]","urban arteries [2899-5727]"),
                                     include.lowest= TRUE,
                                     right =FALSE)
#### natural breaks:  81, 230, 582
local_connectivity_canoas_post$centrality_fct <- cut(local_connectivity_canoas_post$centrality,
                                      breaks=c(
                                              natural_jenks_breaks_cmc_pre_canoas[1],
                                              natural_jenks_breaks_cmc_pre_canoas[2],
                                              natural_jenks_breaks_cmc_pre_canoas[3],
                                              natural_jenks_breaks_cmc_pre_canoas[4],
                                              natural_jenks_breaks_cmc_pre_canoas[5]),
                                      labels =c("low [1-533]","intermediate [533-1510]","high [1510-2899]","urban arteries [2899-5727]"),
                                      include.lowest= TRUE,
                                      right =FALSE)
### 
library(scales)
local_connectivity_canoas$scaled_lwd <- rescale(local_connectivity_canoas$centrality, to = c(0.5, 8)) 
local_connectivity_canoas_pre_map <- mapview::mapview(local_connectivity_canoas,
                                       zcol="centrality_fct",
                                       lwd = local_connectivity_canoas$scaled_lwd,
                                       col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                       col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                       map.types=c("OpenStreetMap","CartoDB.DarkMatter"),
                                       layer.name ="IC connectivity Pre-Event",
                                       popup=leafpop::popupTable(local_connectivity_canoas, 
                                                        zcol=c("centrality","bidirectid"))) 
local_connectivity_canoas_post_map <- mapview::mapview(local_connectivity_canoas_post,
                                       zcol="centrality_fct",
                                       lwd = local_connectivity_canoas$scaled_lwd,
                                       col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                       col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"),
                                       map.types=c("OpenStreetMap","CartoDB.DarkMatter"),
                                       layer.name ="IC connectivity Post-Event",
                                       popup=leafpop::popupTable(local_connectivity_canoas_post, 
                                                        zcol=c("centrality","bidirectid"))) 

local_connectivity_canoas_pre_map | local_connectivity_canoas_post_map + mapview::mapview(flood_extent,
                                                   col.regions="#dbc3b0",
                                                   map.types=c("OpenStreetMap","CartoDB.DarkMatter"),
                                                   alpha.regions= 0.5,
                                                   layer.name="Flooding layer")

The code chunk [ECR14] shows how to calculate the IntraCity (IC) connectivity combining SQL queries with the R library glue (Hester and Bryan 2024). We used the same central query in line 20 that calculated the least-cost paths. However, we created a routable network for each municipality with the iteration between the line 7 and 40.

[ECR 14] How to calculate edge betweenness centrality in the IntraCity Area

## Table with networks
library(sf)
table_rq_1_2_net <- st_read(connection, "table_rq_1_2_net")
table_snapped <- st_read(connection, "table_snapped")
municipalities_ghs <- st_read(connection, "municipalities_ghs")
  
tic()  
for(idx in 1:9){
  tic(glue("run {idx}/9")) 
  table_name_1_centrality <- glue("net_largest_{as.character(idx)}")
  temp_table_name <- glue("temp_centrality_{idx}")
  table_bidirect_id <- glue("t_bidirect_id_{idx}")
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`temp_table_name`}", .con=connection),
          conn=connection)  
query_1_central <- glue_sql("CREATE TABLE {`temp_table_name`} AS
 SELECT   
 b.id,
 b.geom,
 count(b.geom) AS centrality 
 FROM  pgr_dijkstra('SELECT id,
                            fromid AS source,
                            toid AS target,
                            weight AS cost
                      FROM {`table_name_1_centrality`}',
                      ARRAY(SELECT 
                        snapped_id AS start_id 
                        FROM table_snapped 
                        WHERE shapename_id = {idx}::text),
                      ARRAY(SELECT
                        snapped_id AS end_id 
                        FROM table_snapped
                        WHERE shapename_id = {idx}::text),
                      directed := TRUE) j
                      LEFT JOIN {`table_name_1_centrality`} AS b
                      ON j.edge = b.id
                      GROUP BY  b.id, geom;", .con=connection)

dbExecute(con=connection, query_1_central)
toc()
}
toc()
# Time: 53805 ms

The following code repeated the ECR-14 code for the scenario after the flooding.

How to calculate edge betweenness centrality in the IntraCity Area after the flooding

## Table with networks
library(sf)
table_rq_1_2_net_post <- st_read(connection, "table_rq_1_2_net_post")
table_snapped <- st_read(connection, "table_snapped")
municipalities_ghs <- st_read(connection, "municipalities_ghs")

tic()  
for(idx in 1:9){
  tic(glue("run {idx}/9")) 
  table_name_1_centrality_post <- glue("net_largest_post_{as.character(idx)}")
  temp_table_name_post <- glue("temp_centrality_post_{idx}")
  table_bidirect_id_post <- glue("t_bidirect_id_post_{idx}")
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`temp_table_name_post`}", .con=connection),
            conn=connection)  
  query_1_central <- glue_sql("CREATE TABLE {`temp_table_name_post`} AS
 SELECT   
 b.id,
 b.geom,
 count(b.geom) AS centrality 
 FROM  pgr_dijkstra('SELECT id,
                            fromid AS source,
                            toid AS target,
                            weight AS cost
                      FROM {`table_name_1_centrality_post`}',
                      ARRAY(SELECT 
                        snapped_id AS start_id 
                        FROM table_snapped 
                        WHERE shapename_id = {idx}::text),
                      ARRAY(SELECT
                        snapped_id AS end_id 
                        FROM table_snapped
                        WHERE shapename_id = {idx}::text),
                      directed := TRUE) j
                      LEFT JOIN {`table_name_1_centrality_post`} AS b
                      ON j.edge = b.id
                      GROUP BY  b.id, geom;", .con=connection)
  
  dbExecute(con=connection, query_1_central)
  toc()
}
toc()

For each of the municipalities, we repeated the codes ACR11 to ACR13. The following code created the IC connectivity after the flooding in Canoas (id = 3).

How to determine the final IC connectivity for the municipality of Canoas


CREATE TABLE centrality_bidirect_temp_centrality_post_3 AS
SELECT t1.*,
    t2."bdrctid"
FROM
    temp_centrality_post_3 t1
JOIN
    net_largest_post_3 t2
ON
    t1.id = t2.id;
--- The final product requries  99127   
select max("bdrctid") from centrality_bidirect_temp_centrality_post_3 ;
create sequence bididirect_id_temp_centrality_post_3_seq start 99127;
update centrality_bidirect_temp_centrality_post_3 
set "bdrctid" = nextval('bididirect_id_temp_centrality_post_3_seq')
where "bdrctid" is null ;
---- verify
select count(*) 
from centrality_bidirect_temp_centrality_post_3
where "bdrctid" is null; --- 0
----
create table centrality_bidirect_temp_centrality_post_3_bidirect_group as 
select 
       "bdrctid",
       sum(centrality) as  centrality
from centrality_bidirect_temp_centrality_post_3 
group by "bdrctid"; --- this sum the centrality for duplicated bidirect id
---
---- now recover the id
create table centrality_bidirect_temp_centrality_post_3_bidirect_group_id as 
select t1.*, t2.id
from centrality_bidirect_temp_centrality_post_3_bidirect_group t1
join centrality_bidirect_temp_centrality_post_3 t2 
on t1."bdrctid" = t2."bdrctid"; 
---- add geometries
create table centrality_bidirect_temp_centrality_post_3_cleaned as
select t1.*,
    t2.geom,
    t2.toid,
    t2.fromid,
    t2.weight,
    t2."undrcti"
from 
    centrality_bidirect_temp_centrality_post_3_bidirect_group_id t1
join
    net_largest_post_3 t2
on
  t1.id = t2.id;

SELECT * FROM net_largest_post_3;

6.2 Research question 2

Which municipalities were most affected by the flood based on the edge betweenness centrality metrics?

The Table 6.1 aggregate the centrality values for each municipalities. Nova Santa Rita is the most affected municipality when the CM centrality is assessed. However, this same municipality ranked as 6th in terms of loss of IC connectivity.

How to create an interactive table that show municipalities loss of centrality

library(DT)
 sf::st_drop_geometry(connectivity_municipality_lc_rq) |>
  left_join(sf::st_drop_geometry(connectivity_municipality_cmc_rq),
            by=join_by(gid)) |> 
  mutate(diff_ic = (round((post_lc -pre_lc)/pre_lc*100)),
         diff_cmc = (round((post_gc - pre_gc)/pre_gc*100)),
         pop_aff = paste0(pop_aff.x/1000,"K (", round(pop_aff.x/pop.x*100),"%)"),
         gdp = paste0(pib_mun.x/1000,"K")) |>
  select(municipality.x, pop_aff, diff_cmc, diff_ic,gdp) |> 
DT::datatable(colnames=c("Municipality","Pop.Aff","CMC","IC","GDP(R$"),
              extensions = 'Buttons',
              filter = "top",
              options=list(
                          dom = 'Bfrtip',
                          columnDefs = list(list(className = 'dt-center', targets = c(1:4))),
                          pageLength = 5,
                           buttons = c('copy', 'csv')),
              escape = FALSE)

Table 6.1: Pop.Aff: Population affected. CMC: Core Metropolitan Connectivity. IC: IntraCity connectivity. GDP: Gross domestic product per capita (2021)

The Figure 6.3 contains two layers to visualize the IC and CM connectivity. While Nova Santa Rita was the municipality most affected by the flood with a 99% loss in the CM connectivity, Canoas is the most affected when using the IC connectivity with a 60% loss.

How to visualise a comparison of the centrality before and after the floods by municipalities

library(scales)
library(RColorBrewer)
library(mapview)
# Define a red-to-green palette
palette_red_green <- colorRampPalette(c("red", "yellow", "green"))

connectivity_municipality_cmc_rq$diff_cmc <-  round(((connectivity_municipality_cmc_rq$post_gc - connectivity_municipality_cmc_rq$pre_gc) / connectivity_municipality_cmc_rq$pre_gc*100))
connectivity_municipality_lc_rq$diff_lc <-  round(((connectivity_municipality_lc_rq$post_lc - connectivity_municipality_lc_rq$pre_lc) / connectivity_municipality_lc_rq$pre_lc*100))

municipalities_pre_map <- mapview::mapview(connectivity_municipality_cmc_rq,
                                       zcol="diff_cmc",
                                       map.types=c("OpenStreetMap","CartoDB.DarkMatter"),
                                       layer.name ="Core Metropoplitan Connectivity (CM)",
                                       col.regions = palette_red_green(6),
                                       popup=leafpop::popupTable(connectivity_municipality_cmc_rq, 
                                                        zcol=c("municipality","diff","pop_aff","pib_mun"))) 
municipalities_post_map <-  mapview::mapview(connectivity_municipality_lc_rq,
                                       zcol="diff_lc",
                                       map.types=c("OpenStreetMap","CartoDB.DarkMatter"),
                                       layer.name ="IntraCity Connectivity (IC)",
                                       col.regions = palette_red_green(6),
                                       hide=TRUE,
                                       popup=leafpop::popupTable(connectivity_municipality_lc_rq, 
                                                        zcol=c("municipality","diff","pop_aff","pib_mun"))) 

municipalities_pre_map + (municipalities_post_map + mapview(flood_extent, hide =TRUE, col.regions="#dbc3b0", alpha.regions=.5) )

Figure 6.3: A comparison of the centrality before and after the floods aggregated by municipalities

# Calculate connectivity ## Research question 1 > How is road connectivity affected in pre-disaster and post-disaster scenarios based on edge betweenness centrality? ### Core Metropolitan scale @fig-cmc shows an interactive view of the connectivity before and after the flooding event in the Core Metropolitan area. On the left side, the urban arteries such as BR-290 freeway^[[https://prefeitura.poa.br/eptc/noticias/cerca-de-24-mil-veiculos-passam-pelo-corredor-humanitario-no-primeiro-dia](https://prefeitura.poa.br/eptc/noticias/cerca-de-24-mil-veiculos-passam-pelo-corredor-humanitario-no-primeiro-dia)] is displayed in purple previous to the flood event. On the right side, the flood reduced the connectivity of the urban arteries. Municipalities such as Nova Santa Rita were isolated from the core metropolitan area of Porto Alegre. ```{r} #| echo: false #| eval: true #| warning: false #| message: false #| output: false # Load all the rmarkdown_data ending in gpkg saving you time library(tictoc) library(tidyverse) library(sf) tic() rmarkdown_data <- list.files(path= "/home/ricardors/agile-gscience-2024-rs-flood/data/rmarkdown_data", pattern="\\.gpkg") for (i in seq_along(rmarkdown_data)) { file_name <- tools::file_path_sans_ext(rmarkdown_data[i]) assign(file_name, sf::read_sf(paste0("~/agile-gscience-2024-rs-flood/data/rmarkdown_data/", rmarkdown_data[i]))) } # Load non-gpkg files centrality_both_sankey_v3 <- read.csv("~/agile-gscience-2024-rs-flood/data/rmarkdown_data/centrality_both_sankey_v3.csv") df_centrality_both_sankey <- read.csv("~/agile-gscience-2024-rs-flood/data/rmarkdown_data/df_centrality_both_sankey.csv") # Work on flat sf::sf_use_s2(FALSE) # Load gpkg files for RQ results_data <-list.files(path= "/home/ricardors/agile-gscience-2024-rs-flood/data/results_data", pattern="\\.gpkg") for (i in seq_along(results_data)) { file_name <- tools::file_path_sans_ext(results_data[i]) assign(file_name, sf::read_sf(paste0("~/agile-gscience-2024-rs-flood/data/results_data/", results_data[i]))) } toc() ``` ::: {#fig-cmc} <iframe width="780" height="500" src="data/media/pre_post_cm_centrality.html" title="Quarto Documentation"></iframe> A comparison of the centrality before and after the floods in the core metropolitan (CM) area ::: ```{r} #| eval: false #| echo: true #| code-summary: "How to visualise CM connectivity" library(mapview) library(leaflet) library(leafpop) # color: https://gis.stackexchange.com/questions/312027/r-and-mapview-how-can-i-chose-different-colors-for-different-layers ### For pre-event: 1644, 468, 142 pal <- c("#b3cde2","#8a95c5","#8854a7","#7a0f76") centrality_pre$centrality_fct <- cut(centrality_pre$centrality, breaks=c(0,725,2114,4870,11205), labels =c("low [1-726]","intermediate [726-2114]","high [2114-4870]","urban arteries [4780-1125]"), include.lowest= TRUE, right =FALSE) #### natural breaks: 81, 230, 582 centrality_post$centrality_fct <- cut(centrality_post$centrality, breaks=c(0,725,2114,4870,11205), labels =c("low [1-726]","intermediate [726-2114]","high [2114-4870]","urban arteries [4780-1125]"), include.lowest= TRUE, right =FALSE) ### library(scales) centrality_pre$scaled_lwd <- rescale(centrality_pre$centrality, to = c(0.5, 8)) centrality_pre_map <- mapview::mapview(centrality_pre, zcol="centrality_fct", lwd = centrality_pre$scaled_lwd, col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), map.type=c("OpenStreetMap","CartoDB.DarkMatter"), layer.name ="Centrality Pre-Event", popup=leafpop::popupTable(centrality_pre, zcol=c("centrality","bidirectId"))) centrality_post_map <- mapview::mapview(centrality_post, zcol="centrality_fct", col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), map.type=c("OpenStreetMap","CartoDB.DarkMatter"), lwd = centrality_pre$scaled_lwd, layer.name ="Centrality Post-Event", popup=leafpop::popupTable(centrality_post, zcol=c("centrality","bidirectId"))) centrality_pre_map | centrality_post_map + mapview::mapview(flood_extent, col.regions="#dbc3b0", map.type=c("OpenStreetMap","CartoDB.DarkMatter"), alpha.regions= 0.5, layer.name="Flooding layer") ``` We used pgrouting to calculate the centrality by counting the least-cost paths using pgr_dijkstra between origin and destinations (OD). The code [ECR 10] is based on Daniel i. Patterson post ^[[https://urbandatacyclist.com/2020/04/18/how-to-measure-centrality-among-bike-share-trips-using-pgrouting/](https://urbandatacyclist.com/2020/04/18/how-to-measure-centrality-among-bike-share-trips-using-pgrouting/)] and Matt [@Forrest2023-yy] The line 7 executed the function pgr_dijkstra from pgrouting creating a table of the shortest paths from the 300 OD contained in the table _od_300_snapped_v2_. Later, the line 18 added the spatial geometry from the road network saved in the table _ghs_aoi_net_largest_. A higher number of OD returned the error in handling memory allocation ^[[https://discourse.osgeo.org/t/handling-memory-allocation-with-pgrouting-in-centrality-analysis/110407/4](https://discourse.osgeo.org/t/handling-memory-allocation-with-pgrouting-in-centrality-analysis/110407/4)], so we limited to 300 origins and destinations. However, the visual inspection of the spatial connectivity patterns were very similar between 100 and 300 OD revealing the BR-290 freeway as the main urban arteria in both cases. ```{sql} #| eval: false #| echo: true #| code-summary: "[ACR 10] How to calculate edge betweenness centrality in the Core Metropolitan Area" EXPLAIN ANALYZE CREATE TABLE centrality_ghs_aoi_pre AS SELECT b.id, b.geom, count(geom) as centrality FROM pgr_dijkstra( 'SELECT id AS id, "fromId" AS source, "toId" AS target, weight as cost FROM ghs_aoi_net_largest', ARRAY(SELECT snapped_id AS start_id FROM od_300_snapped_v2), ARRAY(SELECT snapped_id AS end_id FROM od_300_snapped_v2), directed:=TRUE) j LEFT JOIN ghs_aoi_net_largest AS b ON j.edge=b.id GROUP BY b.id, b.geom; --- Time: 80542.531 ms ``` The code chunks from ACR11 to ACR13 handles a characteristic from exporting the ORS road network into R, the bidirectid variable. ```{sql} #| eval: false #| echo: true #| code-summary: "[ACR 11] How to add a bidirectId when not present" CREATE TABLE centrality_ghs_aoi_bidirect_pre AS SELECT t1.*, t2."bdrctId" FROM centrality_ghs_aoi_pre AS t1 JOIN ghs_aoi_net_largest AS t2 ON t1.id = t2.id; --- Time: 164.685 ms select max("bdrctId") from centrality_ghs_aoi_bidirect_pre ; ---99527 create sequence bidirectid_centrality_ghs_aoi start 99527; update centrality_ghs_aoi_bidirect_pre set "bdrctId" = nextval('bidirectid_centrality_ghs_aoi') where "bdrctId" is null ; --- checking that bdrctId is not null select count(*) from centrality_ghs_aoi_bidirect_pre where "bdrctId" is null; ---0 ``` We summed up the centrality values of the roads with same bidirectid. ```{sql} #| eval: false #| echo: true #| code-summary: "[ACR 12] How to use bidirectid to calculcalte centrality" create table centrality_ghs_aoi_bidirect_pre_group as select "bdrctId", sum(centrality) as centrality from centrality_ghs_aoi_bidirect_pre group by "bdrctId"; --- Time: 65.993 ms create table centrality_ghs_aoi_bidirect_pre_group_id as select t1.*, t2.id from centrality_ghs_aoi_bidirect_pre_group t1 join centrality_ghs_aoi_bidirect_pre t2 on t1."bdrctId" = t2."bdrctId"; ``` ```{sql} #| eval: false #| echo: true #| code-summary: "[ACR 13] How to calculate final edge betweenness centrality in the Core Metropolitan Area" create table centrality_ghs_aoi_bidirect_pre_cleaned as select t1.*, t2.geom, t2."toId", t2."fromId", t2.weight, t2."undrctI" from centrality_ghs_aoi_bidirect_pre_group_id t1 join ghs_aoi_net_largest t2 on t1.id = t2.id; --- Time: 565.951 ms ``` This process is repeated again changing the table containing the road network _ghs_aoi_net_largest_ to _ghs_aoi_net_largest_post_. While the table _ghs_aoi_net_largest_ represents the road network of the AoI before the flood, the table _ghs_aoi_net_largest_post_ shows the differences on the road network after the flood. ### Intracity scale @fig-lc shows the changes on the centrality when we limited the AoI to each municipality, in this case the municipality of Canoas. The flood hit the west side of Canoas affecting the Mathias Velho dam ^[[https://prefeitura.poa.br/eptc/noticias/cerca-de-24-mil- veiculos-passam-pelo-corredor-humanitario-no-primeiro-dia](https://prefeitura.poa.br/eptc/noticias/cerca-de-24-mil- veiculos-passam-pelo-corredor-humanitario-no-primeiro-dia)] causing the loss of urban arterias such as Rua Florianópolis. ::: {#fig-lc} <iframe width="780" height="500" src="data/media/pre_post_ic_centrality.html" title="Quarto Documentation"></iframe> A comparison of the centrality before and after the floods in the core metropolitan area ::: ```{r} #| echo: true #| eval: false #| message: false #| warning: false #| code-summary: "How to visualiza a comparison of the centrality before and after the floods in the intracity (IC) area" library(leaflet) library(mapview) library(leafpop) ### For pre-event: 1644, 468, 142 pal <- c("#b3cde2","#8a95c5","#8854a7","#7a0f76") local_connectivity_canoas <- filter(local_connectivity_pre, shapename=="Canoas") local_connectivity_canoas_post <- filter(local_connectivity_post, shapename=="Canoas") library(BAMMtools) natural_jenks_breaks_cmc_pre_canoas <- BAMMtools::getJenksBreaks(local_connectivity_canoas$centrality, k=5) local_connectivity_canoas$centrality_fct <- cut(local_connectivity_canoas$centrality, breaks=c( natural_jenks_breaks_cmc_pre_canoas[1], natural_jenks_breaks_cmc_pre_canoas[2], natural_jenks_breaks_cmc_pre_canoas[3], natural_jenks_breaks_cmc_pre_canoas[4], natural_jenks_breaks_cmc_pre_canoas[5]), labels =c("low [1-533]","intermediate [533-1510]","high [1510-2899]","urban arteries [2899-5727]"), include.lowest= TRUE, right =FALSE) #### natural breaks: 81, 230, 582 local_connectivity_canoas_post$centrality_fct <- cut(local_connectivity_canoas_post$centrality, breaks=c( natural_jenks_breaks_cmc_pre_canoas[1], natural_jenks_breaks_cmc_pre_canoas[2], natural_jenks_breaks_cmc_pre_canoas[3], natural_jenks_breaks_cmc_pre_canoas[4], natural_jenks_breaks_cmc_pre_canoas[5]), labels =c("low [1-533]","intermediate [533-1510]","high [1510-2899]","urban arteries [2899-5727]"), include.lowest= TRUE, right =FALSE) ### library(scales) local_connectivity_canoas$scaled_lwd <- rescale(local_connectivity_canoas$centrality, to = c(0.5, 8)) local_connectivity_canoas_pre_map <- mapview::mapview(local_connectivity_canoas, zcol="centrality_fct", lwd = local_connectivity_canoas$scaled_lwd, col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), map.types=c("OpenStreetMap","CartoDB.DarkMatter"), layer.name ="IC connectivity Pre-Event", popup=leafpop::popupTable(local_connectivity_canoas, zcol=c("centrality","bidirectid"))) local_connectivity_canoas_post_map <- mapview::mapview(local_connectivity_canoas_post, zcol="centrality_fct", lwd = local_connectivity_canoas$scaled_lwd, col.regions=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), col=list("#b3cde2","#8a95c5","#8854a7","#7a0f76"), map.types=c("OpenStreetMap","CartoDB.DarkMatter"), layer.name ="IC connectivity Post-Event", popup=leafpop::popupTable(local_connectivity_canoas_post, zcol=c("centrality","bidirectid"))) local_connectivity_canoas_pre_map | local_connectivity_canoas_post_map + mapview::mapview(flood_extent, col.regions="#dbc3b0", map.types=c("OpenStreetMap","CartoDB.DarkMatter"), alpha.regions= 0.5, layer.name="Flooding layer") ``` The code chunk [ECR14] shows how to calculate the IntraCity (IC) connectivity combining SQL queries with the R library glue [@glue]. We used the same central query in line 20 that calculated the least-cost paths. However, we created a routable network for each municipality with the iteration between the line 7 and 40. ```{r} #| eval: false #| echo: true #| code-summary: "[ECR 14] How to calculate edge betweenness centrality in the IntraCity Area" ## Table with networks library(sf) table_rq_1_2_net <- st_read(connection, "table_rq_1_2_net") table_snapped <- st_read(connection, "table_snapped") municipalities_ghs <- st_read(connection, "municipalities_ghs") tic() for(idx in 1:9){ tic(glue("run {idx}/9")) table_name_1_centrality <- glue("net_largest_{as.character(idx)}") temp_table_name <- glue("temp_centrality_{idx}") table_bidirect_id <- glue("t_bidirect_id_{idx}") dbExecute(glue_sql("DROP TABLE IF EXISTS {`temp_table_name`}", .con=connection), conn=connection) query_1_central <- glue_sql("CREATE TABLE {`temp_table_name`} AS SELECT b.id, b.geom, count(b.geom) AS centrality FROM pgr_dijkstra('SELECT id, fromid AS source, toid AS target, weight AS cost FROM {`table_name_1_centrality`}', ARRAY(SELECT snapped_id AS start_id FROM table_snapped WHERE shapename_id = {idx}::text), ARRAY(SELECT snapped_id AS end_id FROM table_snapped WHERE shapename_id = {idx}::text), directed := TRUE) j LEFT JOIN {`table_name_1_centrality`} AS b ON j.edge = b.id GROUP BY b.id, geom;", .con=connection) dbExecute(con=connection, query_1_central) toc() } toc() # Time: 53805 ms ``` ```{sql} #| eval: false #| echo: false ---- 3 CREATE TABLE centrality_bidirect_temp_centrality_3 AS SELECT t1.*, t2."bidirectid" FROM temp_centrality_3 t1 JOIN net_largest_3 t2 ON t1.fid = t2.fid; --- The final product requries 79941 select max("bidirectid") from centrality_bidirect_temp_centrality_3 ; create sequence bididirect_id_temp_centrality_3_seq start 99494; update centrality_bidirect_temp_centrality_3 set "bidirectid" = nextval('bididirect_id_temp_centrality_3_seq') where "bidirectid" is null ; ---- verify select count(*) from centrality_bidirect_temp_centrality_3 where "bidirectid" is null; --- 0 ---- create table centrality_bidirect_temp_centrality_3_bidirect_group as select "bidirectid", sum(centrality) as centrality from centrality_bidirect_temp_centrality_3 group by "bidirectid"; --- this sum the centrality for duplicated bidirect id --- ---- now recover the id create table centrality_bidirect_temp_centrality_3_bidirect_group_id as select t1.*, t2.fid from centrality_bidirect_temp_centrality_3_bidirect_group t1 join centrality_bidirect_temp_centrality_3 t2 on t1."bidirectid" = t2."bidirectid"; ---- add geometries create table centrality_bidirect_temp_centrality_3_cleamed as select t1.*, t2.the_geom, t2.toid, t2.fromid, t2.weight, t2."unidirectid" from centrality_bidirect_temp_centrality_3_bidirect_group_id t1 join net_largest_3 t2 on t1.fid = t2.fid; ``` The following code repeated the ECR-14 code for the scenario after the flooding. ```{r} #| message: false #| warning: false #| eval: false #| echo: true #| code-summary: "How to calculate edge betweenness centrality in the IntraCity Area after the flooding " ## Table with networks library(sf) table_rq_1_2_net_post <- st_read(connection, "table_rq_1_2_net_post") table_snapped <- st_read(connection, "table_snapped") municipalities_ghs <- st_read(connection, "municipalities_ghs") tic() for(idx in 1:9){ tic(glue("run {idx}/9")) table_name_1_centrality_post <- glue("net_largest_post_{as.character(idx)}") temp_table_name_post <- glue("temp_centrality_post_{idx}") table_bidirect_id_post <- glue("t_bidirect_id_post_{idx}") dbExecute(glue_sql("DROP TABLE IF EXISTS {`temp_table_name_post`}", .con=connection), conn=connection) query_1_central <- glue_sql("CREATE TABLE {`temp_table_name_post`} AS SELECT b.id, b.geom, count(b.geom) AS centrality FROM pgr_dijkstra('SELECT id, fromid AS source, toid AS target, weight AS cost FROM {`table_name_1_centrality_post`}', ARRAY(SELECT snapped_id AS start_id FROM table_snapped WHERE shapename_id = {idx}::text), ARRAY(SELECT snapped_id AS end_id FROM table_snapped WHERE shapename_id = {idx}::text), directed := TRUE) j LEFT JOIN {`table_name_1_centrality_post`} AS b ON j.edge = b.id GROUP BY b.id, geom;", .con=connection) dbExecute(con=connection, query_1_central) toc() } toc() ``` For each of the municipalities, we repeated the codes ACR11 to ACR13. The following code created the IC connectivity after the flooding in Canoas (id = 3). ```{sql} #| eval: false #| echo: true #| code-summary: "How to determine the final IC connectivity for the municipality of Canoas" CREATE TABLE centrality_bidirect_temp_centrality_post_3 AS SELECT t1.*, t2."bdrctid" FROM temp_centrality_post_3 t1 JOIN net_largest_post_3 t2 ON t1.id = t2.id; --- The final product requries 99127 select max("bdrctid") from centrality_bidirect_temp_centrality_post_3 ; create sequence bididirect_id_temp_centrality_post_3_seq start 99127; update centrality_bidirect_temp_centrality_post_3 set "bdrctid" = nextval('bididirect_id_temp_centrality_post_3_seq') where "bdrctid" is null ; ---- verify select count(*) from centrality_bidirect_temp_centrality_post_3 where "bdrctid" is null; --- 0 ---- create table centrality_bidirect_temp_centrality_post_3_bidirect_group as select "bdrctid", sum(centrality) as centrality from centrality_bidirect_temp_centrality_post_3 group by "bdrctid"; --- this sum the centrality for duplicated bidirect id --- ---- now recover the id create table centrality_bidirect_temp_centrality_post_3_bidirect_group_id as select t1.*, t2.id from centrality_bidirect_temp_centrality_post_3_bidirect_group t1 join centrality_bidirect_temp_centrality_post_3 t2 on t1."bdrctid" = t2."bdrctid"; ---- add geometries create table centrality_bidirect_temp_centrality_post_3_cleaned as select t1.*, t2.geom, t2.toid, t2.fromid, t2.weight, t2."undrcti" from centrality_bidirect_temp_centrality_post_3_bidirect_group_id t1 join net_largest_post_3 t2 on t1.id = t2.id; SELECT * FROM net_largest_post_3; ``` ## Research question 2 > Which municipalities were most affected by the flood based on the edge betweenness centrality metrics? The @tbl-municipalities aggregate the centrality values for each municipalities. Nova Santa Rita is the most affected municipality when the CM centrality is assessed. However, this same municipality ranked as 6th in terms of loss of IC connectivity. ```{r} #| code-summary: "How to create an interactive table that show municipalities loss of centrality" #| echo: true #| message: false #| waning: false #| label: tbl-municipalities #| tbl-cap: "Pop.Aff: Population affected. CMC: Core Metropolitan Connectivity. IC: IntraCity connectivity. GDP: Gross domestic product per capita (2021)" library(DT) sf::st_drop_geometry(connectivity_municipality_lc_rq) |> left_join(sf::st_drop_geometry(connectivity_municipality_cmc_rq), by=join_by(gid)) |> mutate(diff_ic = (round((post_lc -pre_lc)/pre_lc*100)), diff_cmc = (round((post_gc - pre_gc)/pre_gc*100)), pop_aff = paste0(pop_aff.x/1000,"K (", round(pop_aff.x/pop.x*100),"%)"), gdp = paste0(pib_mun.x/1000,"K")) |> select(municipality.x, pop_aff, diff_cmc, diff_ic,gdp) |> DT::datatable(colnames=c("Municipality","Pop.Aff","CMC","IC","GDP(R$"), extensions = 'Buttons', filter = "top", options=list( dom = 'Bfrtip', columnDefs = list(list(className = 'dt-center', targets = c(1:4))), pageLength = 5, buttons = c('copy', 'csv')), escape = FALSE) ``` The @fig-mun contains two layers to visualize the IC and CM connectivity. While Nova Santa Rita was the municipality most affected by the flood with a 99% loss in the CM connectivity, Canoas is the most affected when using the IC connectivity with a 60% loss. ```{r} #| echo: true #| warning: false #| message: false #| eval: false #| code-summary: "How to visualise a comparison of the centrality before and after the floods by municipalities" library(scales) library(RColorBrewer) library(mapview) # Define a red-to-green palette palette_red_green <- colorRampPalette(c("red", "yellow", "green")) connectivity_municipality_cmc_rq$diff_cmc <- round(((connectivity_municipality_cmc_rq$post_gc - connectivity_municipality_cmc_rq$pre_gc) / connectivity_municipality_cmc_rq$pre_gc*100)) connectivity_municipality_lc_rq$diff_lc <- round(((connectivity_municipality_lc_rq$post_lc - connectivity_municipality_lc_rq$pre_lc) / connectivity_municipality_lc_rq$pre_lc*100)) municipalities_pre_map <- mapview::mapview(connectivity_municipality_cmc_rq, zcol="diff_cmc", map.types=c("OpenStreetMap","CartoDB.DarkMatter"), layer.name ="Core Metropoplitan Connectivity (CM)", col.regions = palette_red_green(6), popup=leafpop::popupTable(connectivity_municipality_cmc_rq, zcol=c("municipality","diff","pop_aff","pib_mun"))) municipalities_post_map <- mapview::mapview(connectivity_municipality_lc_rq, zcol="diff_lc", map.types=c("OpenStreetMap","CartoDB.DarkMatter"), layer.name ="IntraCity Connectivity (IC)", col.regions = palette_red_green(6), hide=TRUE, popup=leafpop::popupTable(connectivity_municipality_lc_rq, zcol=c("municipality","diff","pop_aff","pib_mun"))) municipalities_pre_map + (municipalities_post_map + mapview(flood_extent, hide =TRUE, col.regions="#dbc3b0", alpha.regions=.5) ) ``` ::: {#fig-mun} <iframe width="780" height="500" src="data/media/pre_post_municipality_centrality.html" title="Quarto Documentation"></iframe> A comparison of the centrality before and after the floods aggregated by municipalities :::