results.qmd

# Results

Key variables: - Centrality --- connectivity - Accessibility ---- vulnerability-resilience

```{r}
#| echo: false
#| output: false
#| warning: false

# Load packages
lapply(c("tidyverse","DT","leaflet","sf","DBI", "RPostgres","dplyr", "mapview","leafpop","leaflet","leafsync","terra","raster","stars", "lwgeom","leaflet.extras2","RColorBrewer","tidygeocoder"),
       require,
       character.only =T)
# Connect to the server
eisenberg_connection <- DBI::dbConnect(RPostgres::Postgres(),
                          user= "docker",
                          password = "docker",
                          host = "localhost",
                          dbname="gis",
                          port = 25432)

## Import data
centrality_weighted_100_bidirect_cleaned <- st_read(eisenberg_connection, "centrality_weighted_100_bidirect_cleaned")
centrality_post <- st_read(eisenberg_connection, "centrality_weighted_100_bidirect_cleaned_post")
centrality_pre <-  st_read(eisenberg_connection, "centrality_weighted_100_bidirect_cleaned")
flooding <- sf::st_read("~/heigit_bookdown/data/flooding_porto.geojson") 
nuts_porto <-st_read(eisenberg_connection, "nuts_porto")
ghs <- st_read(eisenberg_connection, "ghs_aoi_porto_alegre")
hospitals <- st_read(eisenberg_connection, "hospital_rs_destination")
pre_network_largest <- st_read(eisenberg_connection, "porto_alegre_net_largest")
pre_network_component <-st_read(eisenberg_connection, "component_analysis_network_porto") |> sf::st_drop_geometry()
pre_network_component_geom <-st_read(eisenberg_connection, "component_analysis_network_porto") 
pre_network_porto <- st_read(eisenberg_connection, "porto_alegre_net_pre")
centrality_weighted_100_bidirect_cleaned <- st_read(eisenberg_connection,"centrality_weighted_100_bidirect_cleaned") 
centrality_weighted_100_bidirect_cleaned_post <-  st_read(eisenberg_connection, "centrality_weighted_100_bidirect_cleaned_post")
weight_sampling_100_origin_post <-  st_read(eisenberg_connection, "weight_sampling_100_origin_post")
weight_sampling_100_destination_post <-  st_read(eisenberg_connection, "weight_sampling_100_destination_post")
porto_alegre_net_pre_ors <- st_read(eisenberg_connection,"porto_alegre_net_pre_ors")
weight_sampling_100_origin <- st_read(eisenberg_connection, "weight_sampling_100_origin")
dijkstra_resulting_table_pre <- st_read( eisenberg_connection,"dijkstra_resulting_table_pre")
porto_alegre_street_united <- st_read(eisenberg_connection, "porto_alegre_street_united")
flooding_rio_grande_do_sul <- st_read(eisenberg_connection, "flooding_rio_grande_do_sul")
weight_sampling_100_destination <-st_read(eisenberg_connection, "weight_sampling_100_destination") 
prueba_largest_network_post_vertices_pgr <- st_read(eisenberg_connection, "prueba_largest_network_post_vertices_pgr")
porto_alegre_net_largest_vertices_pgr <- st_read(eisenberg_connection, "porto_alegre_net_largest_vertices_pgr")
prueba_largest_network_post <- st_read(eisenberg_connection , "prueba_largest_network_post")
options(scipen=999)
```

## Data preparation

### OSM & ORS data

The area of interest  covered `r round(st_area(st_as_sfc(st_bbox(ghs$geom)))/1000000)` km² of surface including `r nrow(nuts_porto)` municipalities with `r nrow(hospitals)` healthcare facilities. On this area, the road network included `r n_distinct(pre_network_component$component)` components with a total length of `r round(st_length(st_union(pre_network_component_geom$the_geom))/1000)` km. The selected component for the centrality analysis represented the `r round(subset(pre_network_component, component==1, select = distance)/sum(pre_network_component$distance) * 100,2)`% of this total, with a length of `r round(subset(pre_network_component, component==1, select = distance)/1000)` km. The transformed graph network contained `r n_distinct(pre_network_largest$target)` target nodes, `r n_distinct(pre_network_largest$source)` source nodes, and `r n_distinct(pre_network_largest$id)` edges.  Classified as a large settlement [@european_commission_joint_research_centre_ghsl_2023], the estimated population from the GHSl-SMOD was `r round(ghs$pop_2020)` with a total built up area of `r round(ghs$bu_m2_2020/1000000) ` km².  

A total of `r n_distinct(weight_sampling_100_origin_post$pt_id)` origin and `r n_distinct(weight_sampling_100_destination_post$pt_id)` destination generated `r nrow(filter(dijkstra_resulting_table_pre, edge==-1))` shortest paths. The initial naive approach using the regular sampling method generated these points in fields and non-populated areas, which would have skewed the representation of the population distribution. This would have resulted in increasing the centrality of roads that are less relevant to the population's movement patterns during the flooding event. Therefore, using the weighted sampling based on the built-up areas provided a closer model representation of the flooding phenomena. 

After subdividing the flooding extent and adjusting it to the Area of Interest (AoI), the flooding covered `r round(st_area(flooding$geometry)/ (10**3)**2) ` km² instead of the original `r round(st_area(flooding_rio_grande_do_sul$geom)/ (10**3)**2)` km². This operation allowed to create the post-flooding road network with less computing costs by only selecting the flooding extent that intersected with the AoI. This post-flooding network with a length of `r round(st_length(st_union(centrality_pre$the_geom))/1000)` km  included the components 21, 14 and 5187, summing up a total of `r  n_distinct(prueba_largest_network_post_vertices_pgr$id)` vertices against the `r n_distinct(porto_alegre_net_largest_vertices_pgr$id)` vertices from the pre-flooding network. When compared to the network before the flooding, the post-network network reduced its length by `r round((st_length(st_union(prueba_largest_network_post$the_geom)) -  st_length(st_union(pre_network_largest$the_geom))) / st_length(st_union(pre_network_largest$the_geom))*100)` % , representing a loss of `r round((st_length(st_union(prueba_largest_network_post$the_geom)) -  st_length(st_union(pre_network_largest$the_geom)))/1000) ` km.

```{r}
#| eval: false
#| code-summary: "PostGIS: Result from the responsive code to import OSM"
# Osmium command
osmium extract -b -51.2791,-30.1722,-50.9407,-29.8048 sul-240501.osm.pbf -o puerto_alegre_urban_center.osm.pbf
```

```{=html}

<iframe width="760" height="500" src="/media/map_sampling_example_v0.html" title = "Inspect of the components in the graph network "></iframe>
```

```{r}
#| message: false
#| warning: false
#| code-summary: "R: Visualizing regular and weighted sampling"
#| eval: false

## Crop and Reproject
## gdalwarp -te -4850853.201784615 -3737074.296348413 -4616291.881935796 -3495378.804761388 GHS_BUILT_V_E2020_GLOBE_R2023A_54009_100_V1_0.tif GHS_BUILT_V_E2020_GLOBE_R2023A_54009_100_V
## gdalwarp -t_srs "EPSG:4326" GHS_BUILT_V_E2020_GLOBE_R2023A_54009_100_V1_0_RioGrandeDoSul.tif GHS_BUILT_V_E2020_GLOBE_R2023A_4326_100_V1_0_RioGrandeDoSul.tif
pal <- mapview::mapviewPalette("mapviewTopoColors")
ghs_build <- stack("/home/ricardo/heigit_bookdown/data/GHS_BUILT_V_E2020_GLOBE_R2023A_4326_100_V1_0_RioGrandeDoSul.tif")
ghs_smod <- stack("/home/ricardo/heigit_bookdown/data/GHS_SMOD_E2020_GLOBE_R2023A_4326_1000_V2_0_RioGrandeDoSul.tif")
ghs_smod_terra <- terra::rast("/home/ricardo/heigit_bookdown/data/GHS_SMOD_E2020_GLOBE_R2023A_4326_1000_V2_0_RioGrandeDoSul.tif")
regular_sampling <- st_read("/home/ricardo/heigit_bookdown/data/random_points_snapped.geojson")
regular_sampling_100 <- sample_n(regular_sampling, size = 100)
weighted_sampling_origin <- st_read(eisenberg_connection,"weight_sampling_100_origin")
weighted_sampling_destination <- st_read(eisenberg_connection,"weight_sampling_100_destination")
nuts <- st_read(eisenberg_connection, "nuts")

weighted_sampling_origin$sample <- "origin"
weighted_sampling_destination$sample <- "destination"
weighted_sampling_both <- rbind(weighted_sampling_origin,
                                weighted_sampling_destination )

weighted_sampling <- weighted_sampling_both |> mutate(sample = as.factor(weighted_sampling_both$sample))

poi_hospital <- st_read(eisenberg_connection,"hospital_rs_node_v2")
m <- matrix(c(
  0, 10, NA,   # Values >= 0 and <= 10 become 0
  10, 13, 1,  # Values > 10 and <= 21 become 1
  13, 29, 2,  # Values > 21 and <= 29 become 2
  30, 30, 3   # Values == 30 become 3
), ncol = 3, byrow = TRUE)
## Classify using the correct matrix
rc2 <- classify(ghs_smod_terra, m, include.lowest=TRUE)
rc2_factor <- as.factor(rc2)
levels(rc2_factor) <- data.frame(
  ID = 1:3,    # These should match the values in the classification
  category = c("Rural: 10-13", "Suburban: 13-29", "Urban Center: 30")
)
category_colors <- c("#008f44","#dedb96", "#cc9152")
mapview(ghs_build[[1]],
        layer.name ="Built-up volume",
        col.regions = pal(100),
        alpha.regions= 0.45,
        hide=TRUE) +
  mapview(ghs_smod[[1]],
          layer.name = "Settlement classification",
          col.regions = pal(100),
          alpha.regions= 0.35,
          hide=TRUE) +
  mapview(weighted_sampling,
          layer.name="Weighted samples",
          zcol="sample",
          col.regions=c("#2D5CA4","#00A3A0"),
          hide=FALSE,
          cex= 3) +
  mapview(regular_sampling_100,
          color = "darkgray",
          col.regions="darkgray",
          cex= 3,
          legend= FALSE,
          hide=TRUE) +
  mapview(subset(poi_hospital,
                 select=c("cd_cnes",
                          "ds_cnes",
                          "id",
                          "geom_hospital")),
          layer.name = "POI - Hospitals",
          color= "darkred",
          col.regions="#CA2334",
          popup=popupTable(poi_hospital, zcol=c("cd_cnes","ds_cnes","id"))) +
mapview(flooding,
            color="darkblue",
            alpha.regions= 0.5,
        hide =TRUE,
            layer.name="Flooding layer") +
    mapview::mapview(centrality_post,
                     lwd = 0.2,
                     color="#cb2a32",
                     hide = TRUE,
                     layer.name ="Post-flooding centrality network") +
    mapview::mapview(centrality_pre,
                     color = "#00a4a4",
                     lwd= 0.2,
                     hide=TRUE,
                     layer.name ="Pre-flooding centrality network") 

```

```{sql}
#| connection: eisenberg_connection
#| eval: false
#| warning: false
#| message: false
#| echo: false

---- Generate the subset to be visualize

CREATE TABLE porto_alegre_net_largest_subset AS                  
SELECT
  st_intersection(net.the_geom, bbox.geom) 
FROM 
  porto_alegre_net_largest AS net
WHERE 
porto_alegre_net_largest.the_geom && subset_post_scenario_bbox; 
---- Generate the subset

CREATE TABLE subset_post_scenario_bbox AS
SELECT st_setsrid(
            st_makeenvelope(
              -51.214287,-30.020226,-51.12934,-29.945862),4326) AS geom;

--- pre_net_subset
CREATE TABLE porto_alegre_net_largest_subset AS                  
SELECT
  *
FROM 
  porto_alegre_net_largest AS net
WHERE 
net.the_geom && 
st_setsrid(
            st_makeenvelope(
              -51.214287,-30.020226,-51.12934,-29.945862),4326)
--- subset_network_in

------------- network inside

CREATE TEMPORARY TABLE flooding_subdivided_porto_join AS
SELECT st_union(the_geom) AS the_geom FROM flooding_subdivided_porto ;

----
CREATE TABLE porto_alegre_street_in_v3 AS
SELECT net.id,
    CASE 
        WHEN ST_Contains(flood.the_geom, net.the_geom)
        THEN net.the_geom
        ELSE st_intersection(net.the_geom, flood.the_geom)
    END AS  geom
FROM porto_alegre_net_largest_subset net
JOIN flooding_subdivided_porto_join flood
ON st_intersects(net.the_geom, flood.the_geom);         

--- subset_network_out
CREATE TABLE porto_alegre_street_out_v3 AS
SELECT net.*
FROM porto_alegre_net_largest_subset
WHERE net.id NOT IN (
    SELECT net.id
    FROM porto_alegre_street_in_v3 net);
    
---- network outside

CREATE TABLE porto_alegre_net_outside_v3 AS
WITH porto_alegre_ghs AS(
    SELECT 
       ghs.*
    FROM
       urban_center_4326 AS ghs
    JOIN
       nuts 
    ON 
       st_intersects(nuts.geom, ghs.geom)
    WHERE 
        nuts.shapename = 'Porto Alegre'),
--- Bounding Box that contained the GHS in Porto Alegre
porto_alegre_ghs_bbox AS(
    SELECT 
        st_setsrid(st_extent(geom),4326) as geom_bbox
    FROM 
        porto_alegre_ghs),
flooding_sul_subdivided AS (
        SELECT 
            st_subdivide(geom) as the_geom
        FROM
            flooding_rio_grande_do_sul),
exterior_ring_porto_alegre_v2 AS (
SELECT 
    ST_ExteriorRing((ST_Dump(union_geom)).geom) as geom
FROM (
    SELECT 
        ST_Union(flood.the_geom) as union_geom
    FROM 
        porto_alegre_ghs_bbox as bbox
    JOIN 
        flooding_sul_subdivided as flood
    ON 
        ST_Intersects(flood.the_geom, bbox.geom_bbox)
) AS subquery)
SELECT net.id,
    CASE
        WHEN NOT ST_Contains(flood.geom, net.the_geom)
        THEN net.the_geom
            ELSE st_intersection(net.the_geom, flood.geom)
    END AS  geom,
    net.target,
       net.source,
       cost,
       "unidirectid",
       "bidirectid"
FROM
porto_alegre_net_largest_subset AS net
JOIN exterior_ring_porto_alegre_v2 flood ON
st_intersects(net.the_geom, flood.geom)
WHERE 
  net.the_geom && 
    st_setsrid(
    st_makeenvelope(-51.214287,-30.020226,-51.12934,-29.945862),4326);
----
--- For the network
CREATE INDEX idx_porto_alegre_net_outside_v2 ON porto_alegre_net_outside_v2 USING gist (geom);

CLUSTER porto_alegre_net_outside_v2 USING idx_porto_alegre_net_outside_v2;

--- For the flooding mask
CREATE INDEX flooding_sul_subdivided_idx ON flooding_sul_subdivided USING gist (the_geom);

CLUSTER flooding_sul_subdivided USING flooding_sul_subdivided_idx;

---- Before doing difference
VACUUM(FULL, ANALYZE) porto_alegre_net_outside_v2;
VACUUM(FULL, ANALYZE) flooding_sul_subdivided;
----
CREATE INDEX idx_porto_alegre_net_outside_v3 ON porto_alegre_net_outside_v3 USING gist (geom);

CLUSTER porto_alegre_net_outside_v3 USING idx_porto_alegre_net_outside_v3;

--- For the flooding mask
CREATE INDEX flooding_subdivided_porto_join_idx ON flooding_subdivided_porto_join USING gist (the_geom);

CLUSTER flooding_subdivided_porto_join USING flooding_subdivided_porto_join_idx ;

---- Before doing difference
VACUUM(FULL, ANALYZE) porto_alegre_net_outside_v3;
VACUUM(FULL, ANALYZE) flooding_subdivided_porto_join;

----

CREATE TABLE flooding_symple as 
SELECT st_union(geom) as the_geom FROM flooding_cleaned_porto_union_simple;

CREATE INDEX flooding_symple_idx ON flooding_symple USING gist (the_geom);
CLUSTER flooding_symple USING flooding_symple_idx;

CREATE TABLE difference_outside_flood_v4 AS
SELECT net.id,
        target,
        source,
        cost,
        unidirectid,
        bidirectid,
st_difference(net.geom, flood.the_geom) AS the_geom
FROM porto_alegre_net_outside_v3 AS net,
flooding_symple  AS flood;

-----------
CREATE TABLE porto_alegre_street_united_v3 AS
SELECT *
FROM porto_alegre_street_out_v3
UNION
SELECT *
FROM difference_outside_flood_v4;

--- Final product: porto_alegre_street_united_v2
```

## RQ1: Centrality analysis

> How does road connectivity change after being impacted by flooding based on connectivity metrics?

### Sum Edge Betweenness

The average sum edge betweenness before the flooding was `r round(mean(centrality_weighted_100_bidirect_cleaned$centrality))`, being decreased to `r round(mean(centrality_weighted_100_bidirect_cleaned_post$centrality))` after the flooding. In contrast, the flooding caused a `r round((round(max(centrality_weighted_100_bidirect_cleaned_post$centrality)) - round(max(centrality_weighted_100_bidirect_cleaned$centrality))) / round(max(centrality_weighted_100_bidirect_cleaned$centrality)) * 100) ` % decrease on the maximum values,  reducing the sum edge betweenness from `r max(centrality_weighted_100_bidirect_cleaned$centrality)` to `r max(centrality_weighted_100_bidirect_cleaned_post$centrality)`.  Prior to the event, the `r round(round(st_length(st_union((subset(filter(centrality_weighted_100_bidirect_cleaned, centrality > quantile(centrality_weighted_100_bidirect_cleaned$centrality)[4]), select=the_geom)))))/round(st_length(st_union(subset(centrality_weighted_100_bidirect_cleaned, select=the_geom))))*100)`% of the network that extended over `r round(st_length(st_union((subset(filter(centrality_weighted_100_bidirect_cleaned, centrality > quantile(centrality_weighted_100_bidirect_cleaned$centrality)[4]), select=the_geom))))/1000)` km, accounted for `r round(sum(centrality_pre$centrality[centrality_pre$centrality > quantile(centrality_pre$centrality)[4]])/sum(centrality_pre$centrality)*100)` % of the overall connectivity adding `r sum(centrality_pre$centrality[centrality_pre$centrality > quantile(centrality_pre$centrality)[4]])` to the network.  After the event, the upper quartile of the network, with a length of `r round(st_length(st_union((subset(filter(centrality_weighted_100_bidirect_cleaned_post, centrality > quantile(centrality_weighted_100_bidirect_cleaned_post$centrality)[4]), select=the_geom))))/1000)`km, maintained `r round(sum(centrality_post$centrality[centrality_post$centrality > quantile(centrality_post$centrality)[4]])/sum(centrality_post$centrality)*100) `% of the total centrality, however, it only contributed `r  sum(centrality_post$centrality[centrality_post$centrality > quantile(centrality_post$centrality)[4]]) ` to the sum edge betweenness.


```{r}
#| eval: false
#| code-summary: "R: Visualizing changes of centrality on the network before and after the flooding calculated with pgrouting"

### natural breaks
#### For pre-event: 1644, 468, 142
centrality_pre$centrality_fct <- cut(centrality_pre$centrality,
                  breaks=c(0,142,468,1644),
                  labels =c("low","medium","high"),
                  include.lowest= TRUE,
                  right =FALSE)
#### natural breaks:  81, 230, 582
centrality_post$centrality_fct <- cut(centrality_post$centrality,
                  breaks=c(0,81,230,582),
                  labels =c("low","medium","high"),
                  include.lowest= TRUE,
                  right =FALSE)
### 
centrality_pre_map <- mapview::mapview(centrality_pre,
                                       zcol="centrality_fct",
                                        lwd ="centrality",
                                       layer.name ="Centrality Pre-Event",
              popup=popupTable(centrality_pre, 
                               zcol=c("id","centrality","bidirectid"))) 
centrality_post_map <- mapview::mapview(centrality_post,
                                       zcol="centrality_fct",
                                        lwd = "centrality",
                                       layer.name ="Centrality Post-Event",
              popup=popupTable(centrality_pre, 
                               zcol=c("id","centrality","bidirectid"))) 

centrality_pre_map | centrality_post_map + mapview(flooding,
          color="darkblue",
          alpha.regions= 0.5,
          layer.name="Flooding layer")

```

```{=html}

<iframe width="760" height="500" src="/media/centrality_map_results.html" title = "Sum Edge Betweenness in the Urban Settlement located in Porto Alegre "></iframe>
```

The most central road segments were georeferenced to reveal the name of these identified critical roads. It is found that urban arteries such as Freeway Anchieta or Avenida Presidente Castello Branco lost their entire centrality after the flooding. Other roads such as Avenida Teresópolis or Avenida Manoel Elias, Jadim Leopoldina, even being affected by the flooding, only lost their centrality partially.In contrast, roads which were not used before the flooding, such as those located nearby Costa do Morro, Itacolomi or Costa do Ipiranga, were used after the flooding event. These results are shown in the table 1 and mapped in the figure 1.

```{r}
#| message: false
#| warning: false
#| code-summary: "R: Table of the most central road segments georeferenced"
#| eval: false

library(plyr)
## Categories for pre-flooding or post-flooding
df_centrality_pre <- centrality_pre |>
                      sf::st_drop_geometry() |>
                      mutate(event = "pre-flooding") 
df_centrality_post <- centrality_post |>
                        sf::st_drop_geometry() |>
                    mutate(event = "post-flooding")

#### For pre-event: 1644, 468, 142
df_centrality_pre$centrality_fct <- cut(df_centrality_pre$centrality,
                  breaks=c(0,142,468,1644),
                  labels =c("low","medium","high"),
                  include.lowest= TRUE,
                  right =FALSE)
#### natural perk:  81, 230, 582
df_centrality_post$centrality_fct <- cut(df_centrality_post$centrality,
                  breaks=c(0,81,230,582),
                  labels =c("low","medium","high"),
                  include.lowest= TRUE,
                  right =FALSE)
## Join both in one dataframe
df_centrality_both <- left_join(df_centrality_pre,
                    df_centrality_post,
                    by = 'id',
                    suffix=c("_pre","_post")) |>    dplyr::select(c("id","bidirectid_pre","bidirectid_post","centrality_pre","centrality_post","centrality_fct_pre","centrality_fct_post","event_post","event_pre")) |>
  mutate(centrality_post = replace_na(centrality_post, 0)) ## Roads covered by the flooding appeared as NA in the post-scenario, replace tha NA value for 0

## Calculate change on centrality after being impacted by flooding 
df_centrality_both <-  df_centrality_both |>
            mutate(centrality_diff = centrality_post-centrality_pre)
### Obtain distribution using quantiles
df_centrality_both_quantiles <- quantile(df_centrality_both$centrality_diff, na.rm =TRUE)
### Filter the worst scenario, most negative values below second quartile
df_centrality_negative_outlier <- df_centrality_both |>
                        filter( centrality_diff <= df_centrality_both_quantiles[2])
### Join with original data to obtain again geometry
df_centrality_negative_outlier_geom <- df_centrality_negative_outlier |>
                        arrange(centrality_diff) |>
                        slice(1:1000) |>
                        left_join(centrality_pre, by = "id") |>
                    mutate(the_geom_centroid = st_centroid(the_geom),
                           lon = st_coordinates(the_geom_centroid)[,1],
                           lat = st_coordinates(the_geom_centroid)[,2])
# Obtain a subset of the 500 observations oredered by lowest values on centrality difference
df_centrality_negative_outlier_geom_unique <- df_centrality_negative_outlier_geom |> 
  distinct(centrality_pre, .keep_all = TRUE) 
# Reverse geocoding to obtain the address based on the centroids
df_centrality_negative_outlier_distinct_address <-  tidygeocoder::reverse_geocode(df_centrality_negative_outlier_geom_unique, lat=lat, lon=lon, method="osm")
## Create ID for the found addresses (97)
df_centrality_negative_outlier_distinct_address$address_number <- seq(1, 167,1)

df_centrality_negative_outlier_distinct_address$short_address <- str_extract(df_centrality_negative_outlier_distinct_address$address, "^[^,]+, [^,]+") 
df_centrality_negative_outlier_distinct_address$centrality_diff_perc <- round(((df_centrality_negative_outlier_distinct_address$centrality_post - df_centrality_negative_outlier_distinct_address$centrality_pre) / df_centrality_negative_outlier_distinct_address$centrality_pre *100),2)
  
## DT Table
DT::datatable(subset(df_centrality_negative_outlier_distinct_address, select=c("address_number", "short_address","centrality_pre", "centrality_post","centrality_diff","centrality_diff_perc","id","source","target","the_geom")), 
              colnames= c("ID","Address","Pre-Centrality", "Post-Centrality","Diff-Centrality","Diff-Centrality(%)","id_1","source","target","the_geom"),
              filter="top",
              class='compact', rownames=FALSE, escape=FALSE, caption='Data description',
              extensions=c("Buttons",'RowGroup'),
              options=list(
                  order=list(list(5, 'desc'), list(2,'desc')),  # Sort by the first column (index 5)
                  dom="Bfrtip",
                  columnDefs = list(list(visible=FALSE, targets= c(6,7,8,9))),
                  buttons=c("copy", "csv", "pdf"),
                  initComplete = JS(
                      "function(settings, json) {",
                      "$(this.api().table().header()).css({'background-color': '#d50038', 'color': '#fff'});",
                      "}")
              )
) |> 
      DT::formatStyle("centrality_pre",
     background=DT::styleColorBar(range(df_centrality_negative_outlier_distinct_address$centrality_pre), '#ee8b8b'),
                    backgroundSize='98% 88%',
                    backgroundRepeat='no-repeat',
                    backgroundPosition='center') 
```


```{r}
#| eval: false
df_centrality_bind <- bind_rows(df_centrality_pre, df_centrality_post)

df_centrality_bind <- df_centrality_bind |> 
  mutate(centrality_fct = case_when(
    centrality >= 1 & centrality < 53 ~ "[1-53]",  # Change > to >=
    centrality >= 53 & centrality < 117 ~ "[53-117]",
    centrality >= 117 & centrality < 700 ~ "[117-700]",
    centrality >= 700 & centrality < 1644 ~ "[700-1644]",
    centrality == 0 ~ "[0]",
    TRUE ~ NA_character_  
  ))

df_centrality_pre

library(metan)
library(ggsankey)
library(highcharter)

d <- select(df_centrality_bind, c(event, centrality_fct))


hchart(data_to_sankey(d), "sankey", name = "Centrality pre and post flooding") |> 
hc_title(text= "Sankey Diagram: Porto Alegre") %>%
  hc_subtitle(text= "Centrality analysis")  %>%
  hc_caption(text = "<b>Change on the distribution edge betweenness<b>")%>%
  hc_add_theme(hc_theme_economist())

#### secon sankey diagram
library("ggalluvial")
pre_network_largest
sankey_pre_centrality <- centrality_pre |> dplyr::select(c("centrality", "the_geom","event")) |>
                                            mutate(length = st_length(the_geom))
sankey_post_centrality <- centrality_post |> dplyr::select(c("centrality", "the_geom","event")) |> 
                                              mutate(length = st_length(the_geom))
sankey_na_centrality <- pre_network_largest |> mutate(centrality= rep(0, nrow(pre_network_largest)),
                                               event = rep ("not-used", nrow(pre_network_largest)),
                                               length = st_length(the_geom)) |> 
                              dplyr::select(c("centrality", "the_geom","event","length"))
sankey_all_centrality <- rbind(sankey_pre_centrality, sankey_post_centrality,sankey_na_centrality)

sankey_all_centrality$centrality_fct <- cut(sankey_all_centrality$centrality,
                                  breaks = c(1,13,53,117,1644),
                                  labels= c("[1-13]","[13-53]","[53-117]","[117-1644]"),
                                  include.lowest =TRUE)

# Summarize the total length for each event and centrality quartile
sankey_post_pre <-rbind(sankey_pre_centrality, sankey_post_centrality) 
sankey_post_pre$centrality_fct <- cut(sankey_post_pre$centrality,
                                  breaks = c(1,13,53,117,700,2000),
                                  labels= c("[1-13]","[13-53]","[53-117]","[117-300]","[300-2000]"),
                                  include.lowest =TRUE)

sankey_agg <- sankey_post_pre %>% sf::st_drop_geometry() |> 
  dplyr::group_by(event, centrality_fct) %>%
  dplyr::summarise(total_length = sum(as.numeric(length))) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(
    total_length_all = sum(total_length),
    percent = (total_length / total_length_all) * 100
  )

(unlist(dplyr::summarise(group_by(df_sankey_post_pre, event), t_length = sum(length))[1,2])/ unlist(dplyr::summarise(group_by(df_sankey_post_pre, event), t_length = sum(length))[2,2])) * unlist(dplyr::summarise(group_by(df_sankey_post_pre, event), t_length = sum(length))[1,2])
#### oooootro intentooooo 
library(forcats)
library(data.table)
sankey_wider <- sankey_post_pre |> pivot_wider(names_from= event, values_from = centrality)
df_sankey_wider <- sankey_wider |> sf::st_drop_geometry() 
  
df_sankey_wider <- df_sankey_wider |> mutate(`pre-flooding` = 
                            case_when(is.na(df_sankey_wider$`pre-flooding`) ~ 0,
                                       .default=df_sankey_wider$`pre-flooding`),
                          `post-flooding`= case_when(is.na(df_sankey_wider$`post-flooding`) ~ 0,
                                      .default=df_sankey_wider$`post-flooding`)
                          )
df_sankey_wider_zero <- df_sankey_wider |> 
    mutate(pre_flooding_fct = cut(`pre-flooding`,
                                  breaks = c(-Inf, 0, 1, 13, 53, 117, 700, 1645, Inf),
                                  labels = c("[0]", "[1]", "[1-13]", "[13-53]", "[53-117]", "[117-300]", "[300-1645]", "[>1645]"),
                                  include.lowest = TRUE),
           pre_flooding_fct = forcats::fct_explicit_na(pre_flooding_fct, "[0]"),
           post_flooding_fct = cut(`post-flooding`,
                                   breaks = c(-Inf, 0, 1, 13, 53, 117, 700, 1645, Inf),
                                   labels = c("[0]", "[1]", "[1-13]", "[13-53]", "[53-117]", "[117-300]", "[300-1645]", "[>1645]"),
                                   include.lowest = TRUE),
           post_flooding_fct = forcats::fct_explicit_na(post_flooding_fct, "[0]"))



hchart(data_to_sankey(df_sankey_post_pre_wider_v4), "sankey", name = "Centrality pre and post flooding")

df_sankey_wider_zero |> select(c(pre_flooding_fct, post_flooding_fct)) |>  group_by(pre_flooding_fct, post_flooding_fct) |> summary(n=n(), .groups= "drop")

ggplot(data = df_sankey_wider_zero,
       aes(axis1 = pre_flooding_fct, axis2 = post_flooding_fct, 
           y = Freq)) +
  scale_x_discrete(limits = c("Class", "Sex", "Age"), expand = c(.2, .05)) +
  xlab("Demographic") +
  geom_alluvium(aes(fill = Survived)) +
  geom_stratum() +
  geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
  theme_minimal()

library(networkD3)
df_network <- sankey_post_pre |> sf::st_drop_geometry() |> select(c("event","centrality_fct"))
df_network_v2 <- df_network |> group_by(centrality_fct, event) |> summarise(n=n())
df_network_v2 <- df_network_v2 %>%
  bind_rows(tibble(centrality_fct = "[300-2000]",
                   event = "post-flooding",
                   n = 0))

links <- data.frame(source=c(df_network_v2$event),
                    target = c(df_network_v2$centrality_fct),
                    value = c(df_network_v2$n))

nodes <- data.frame(
  name=c(as.character(links$source), 
  as.character(links$target)) %>% unique()
)

links$IDsource <- match(links$source, nodes$name)-1 
links$IDtarget <- match(links$target, nodes$name)-1

## nice gepeto

  

df_long <- df_network_v2 |> 
    ungroup() |> 
    add_row(centrality_fct = "[300-2000]", event = "post-flooding", n = 0) |> 
    pivot_wider(names_from = event, values_from = n) |> 
    mutate(alluvium = row_number())

ggplot(data = df_long,
       aes(axis1 = centrality_fct, axis2 = event, y = n)) +
  geom_alluvium(aes(fill = event)) +
  geom_stratum() +
 geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("Pre-Centrality", "Post-Centrality")) +
  theme_minimal() +
  ggtitle("Centrality Analysis Before and After Flooding",
          "Comparing centrality distributions across events")


ggplot(df_long,
       aes(axis1 = `pre-flooding`, axis2 = alluvium, y = `post-flooding`)) +
  geom_alluvium(aes(fill = centrality_fct)) +
  geom_stratum(width = 1/6) +
  geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("Centrality", "Alluvium"),
                   expand = c(.1, .1)) +
  labs(title = "Alluvial Plot of Centrality Pre- and Post-Flooding",
       y = "Number of Roads", x = "Centrality Ranges") +
  theme_minimal()

df_long <- df_long %>%
  pivot_longer(cols = c(`post-flooding`, `pre-flooding`), 
               names_to = "event", 
               values_to = "n")

df_long <- df_sankey_wider_zero |> 
  pivot_longer(cols = c(`pre-flooding`, `post-flooding`), 
               names_to = "event", 
               values_to = "n")
ggplot(df_long, aes(axis1 = pre_flooding_fct, 
                    axis2 = post_flooding_fct, 
                    y = n)) +
  geom_alluvium(aes(fill = event), width = 1/12) +
  geom_stratum() +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), 
            size = 3, color = "white") +
  scale_x_discrete(limits = c("Pre-Flooding Centrality", "Post-Flooding Centrality"),
                   expand = c(0.05, 0.05)) +
  labs(title = "Alluvial Diagram of Centrality Changes",
       x = "Centrality Categories",
       y = "Number of Roads (n)") +
  theme_minimal() +
  theme(legend.position = "top")

```


```{=html}

<iframe width="760" height="500" src="/media/centrality_table_results.html" title = "Sum Edge Betweenness in the Urban Settlement located in Porto Alegre "></iframe>
```

Regarding the distribution of the centrality, the histogram on the figure 1 represented the betweenness centrality on a log10 scale on the X-axis, while the Y axis shown the frequency of occurence for both scenarios. In both scenarios, the centrality distribution exibited strong positive skewness, with most of the values concentrated at the lower values of the edge betweenness. It is observed that the frequency of mid-range centrality values between 10 and 100 increased after the flooding event, however, the tail with larger values after 1000 was no longer present. The extreme positive centrality values found only in the pre-disaster network affected the general centrality lowering the mean.  

```{r}
#| warming: false
#| message: false
#| code-summary: "R: Histogram of betweenness centrality calculated with pgrouting"

library(plyr)
# Tidy data and wrangling
centrality_post$event <- "post-flooding"
centrality_pre$event <- "pre-flooding"  
centrality_both <- rbind(centrality_post[,c("id","centrality","event")] ,
                                   centrality_pre[,c("id","centrality","event")])
mu <- plyr::ddply(centrality_both, "event", summarise, grp.mean=mean(centrality))

# Create histogram
ggplot(centrality_both, aes(x=centrality, fill=event)) +
                    geom_histogram(alpha=0.4,) + scale_x_log10() +
                  labs(title = "Centrality analysis",
                      subtitle= "Histogram on log10 scale",
                      x = " Betweenness centrality (centrality)",
                      y = "Frequency (count)") +
  theme(plot.title=element_text(family ="bold", hjust=0.5),
        plot.subtitle = element_text(colour="#626262", hjust=0.5),
        legend.position='bottom') +
    theme_minimal()
```

The following grouped bar plot illustrates the sum of the edge betweenness values classified by their quantiles for both events, before and after the flooding. The X-Axis is categorized in the four quantiles, while yhe y-axis summed up the betweenness values of these quartiles. In Both events a similar pattern is found, the roads belonging to the Q4 caused most of the edge beetweenness of the network. The `r round(st_length(st_union(filter(centrality_pre, centrality < quantile(centrality_pre$centrality)[4])["the_geom"]))/1000)` km that comprised all the roads below the fourth quartile only contributed to `r round(summarise(filter(sf::st_drop_geometry(centrality_pre), centrality < quantile(centrality_pre$centrality)[4])["centrality"], centrality = sum(centrality))/sum(sf::st_drop_geometry(centrality_pre)$centrality) * 100)` for the preflooding scenario. Therefore, the findings shown an positive skewed distribution, where a small proportion of roads were responsible for the  roads held responsible for the `r round(summarise(filter(sf::st_drop_geometry(centrality_pre), centrality > quantile(centrality_pre$centrality)[4])["centrality"], centrality = sum(centrality))/sum(sf::st_drop_geometry(centrality_pre)$centrality) * 100)`% of the centrality of the network.


```{r}
#| warming: false
#| message: false
#| code-summary: "R: Betweenness centrality values grouped by their quantiles calculated with pgrouting "

df_centrality_pre <- centrality_pre |>
                        sf::st_drop_geometry() |>
                        mutate(event = "pre-flooding") 
df_centrality_post <- centrality_post |> 
                          sf::st_drop_geometry() |>
                              mutate(event = "post-flooding")
##
df_centrality_both_barplot <- dplyr::bind_rows(df_centrality_pre,
                                               df_centrality_post)
### percentiles
q <- quantile(df_centrality_both_barplot$centrality)
q_pre <- quantile(
  df_centrality_pre[df_centrality_pre$event == 'pre-flooding',]$centrality)
q_post <- quantile(
  df_centrality_post[df_centrality_post$event == 'post-flooding',]$centrality)
### Centrality
df_centrality_both_barplot_cat <- na.omit(df_centrality_both_barplot) |> 
  mutate(centrality_cat=as.factor(case_when(
    centrality <= q_pre[2] ~ "Q1",
    centrality <= q_pre[3] ~ "Q2",
    centrality <= q_pre[4] ~ "Q3",
    centrality >= q_pre[4] ~ "Q4",
    TRUE ~ "missing")))  
## tidy data
barplot_event_cat_sum <- df_centrality_both_barplot_cat |> 
                          dplyr::group_by(event, centrality_cat) |>
                           dplyr::summarise(sum_centrality = sum(centrality))
# basic plot
 barplot_event_cat_sum |>  
  ggplot(aes(x=centrality_cat,
              y=sum_centrality,
              fill=event)) +
  geom_col(width=0.5, position="dodge") +
  labs(title = "Centrality analysis",
       subtitle= "Barplot grouped by the post and pre event",
       x = " Quantiles",
       y = "Edge betweenness") +
  theme(plot.title=element_text(family ="bold", hjust=0.5),
        plot.subtitle = element_text(colour="#626262", hjust=0.5),
        legend.position='bottom') +
     theme_minimal()

```


```{r}
#| warning: false
#| message: false
#| eval: false
#| include: false
#| echo: false
#| code-summary: "R: Identifying"

df_centrality_both_pareto_cat <- na.omit(
  df_centrality_both_barplot) |> 
  mutate(centrality_cat=as.factor(
      case_when(
    centrality <= quantile(
        df_centrality_pre[df_centrality_pre$event == 'pre-flooding',]$centrality, .20) ~ "20%",
    centrality <= quantile(
        df_centrality_pre[df_centrality_pre$event == 'pre-flooding',]$centrality, .40) ~ "40%",
    centrality <= quantile(
        df_centrality_pre[df_centrality_pre$event == 'pre-flooding',]$centrality, .60) ~ "60%",
    centrality <= quantile(
        df_centrality_pre[df_centrality_pre$event == 'pre-flooding',]$centrality, .80) ~ "80%",
    centrality >= quantile(
        df_centrality_pre[df_centrality_pre$event == 'pre-flooding',]$centrality, .80) ~ "80-100%",
    TRUE ~ "missing")))  

d <- df_centrality_both_pareto_cat |>
            st_drop_geometry() |>
            filter(event=="pre-flooding") |>
            dplyr::group_by(centrality_cat) |> dplyr::summarise(centrality_sum=sum(centrality)) |> 
  arrange(desc(centrality_sum)) |>
  mutate(cumsum=cumsum(centrality_sum),
         freq=round(centrality_sum/sum(centrality_sum),3),
         cum_freq=cumsum(freq))
## Saving Parameters 

def_par <- par() 

# New margins
par(mar=c(5,5,4,5)) 

## plot bars, pc will hold x values for bars
pc = barplot(d$centrality_sum,
             width = 1, space = 0.2, border = NA, axes = F,
             ylim = c(0, 1.05 * max(d$centrality_sum, na.rm = T)), 
             ylab = "Frequency of betweenness centrality" ,
             xlab = "Quantiles",
             cex.names = 0.7, 
             names.arg = d$centrality_cat,
             main = "Identifying significant root causes for centrality")

## anotate left axis
axis(side = 2, at = c(0, d$centrality_sum), las = 1, col.axis = "grey62", col = "grey62", tick = T, cex.axis = 0.8)

## frame plot
box( col = "grey62")

## Cumulative Frequency Lines 
px <- d$cum_freq * max(d$centrality_sum, na.rm = T)
lines(pc, px, type = "b", cex = 0.7, pch = 19, col="#d50038")

## Annotate Right Axis
axis(side = 4, at = c(0, px), labels = paste(c(0, round(d$cum_freq * 100)) ,"%",sep=""), 
     las = 1, col.axis = "grey62", col = "#d50038", cex.axis = 0.8, col.axis = "#d50038")

## restoring default paramenter
par(def_par) 

#### ![](/media/pareto_chart.png)
```

## RQ2: Accessibility analysis

> Which healthcare facilities will be most affected by flooding based on accessibility metrics?

### Sum Edge Betweenness:

```{r}
#| eval: false
#| code-summary: "R: Visualizing the accesibility most affected from healtchare facilities based on betweenness centrality"

## Eisenberg
hospitals_betweenness_pre <- st_read(eisenberg_connection, "centrality_424_hospitals_porto_end_id_centrality") 
hospitals_betweenness_post <- st_read(eisenberg_connection,"centrality_424_hospitals_porto_end_id_centrality_post")
### Locally
hospitals_betweenness_pre_local <- st_read(eisenberg_connection, Id(schema="heigit", table = "centrality_424_hospitals_porto_end_id_centrality_group"))
hospitals_betweenness_post_local <- st_read(eisenberg_connection, Id(schema="heigit", table = "centrality_424_hospitals_porto_end_id_centrality_post_group"))
hospitals_betweenness_pre_local$event <- "pre-event"
hospitals_betweenness_post_local$event <- "post-event"
## Import
hospitals_betweenness_both_long <- hospitals_betweenness_pre_local |> bind_rows(hospitals_betweenness_post_local) |> subset(select=c("cd_cnes","ds_cnes","max_centrality","event")) |> pivot_wider(
    names_from = event,
    values_from = max_centrality
  ) |>  rename( betwenness_pre= `pre-event`,
                betwenness_post= `post-event`) 
hospitals_betweenness_both_long[is.na(hospitals_betweenness_both_long)] <- 0    
## pivot
hospitals_betweenness_long <- hospitals_betweenness_both_long |>
                                  mutate(diff_betweenness_percent =
                                           (betwenness_post-betwenness_pre)/(betwenness_pre)*100,
                                         diff = betwenness_post-betwenness_pre,
                                         letter = LETTERS[1:22])  |>
  arrange(letter) |> 
                             pivot_longer(
                                     c(betwenness_pre,betwenness_post),
                                      names_to="betweenness_type",
                                      values_to="betweenness_value") 
                   
hospitals_betweenness_long_pre <- 
  hospitals_betweenness_long |>
            filter(betweenness_type == "betwenness_pre")
hospitals_betweenness_long_post <- 
  hospitals_betweenness_long |>
  filter(betweenness_type == "betwenness_post")
##

df_labels_betweenness <- hospitals_betweenness_long |> 
  group_by(letter) |> 
  summarize(midpoint = mean(betweenness_value), 
            diff_betweenness_percent = first(diff_betweenness_percent))
##
p <- ggplot(hospitals_betweenness_long, aes(x=betweenness_value, y=reorder(letter, desc(diff_betweenness_percent)))) +
  geom_segment(data=hospitals_betweenness_long_pre,
               aes(x=betweenness_value, y = letter,
                   yend=hospitals_betweenness_long_post$letter, xend=hospitals_betweenness_long_post$betweenness_value)) +
  geom_point(aes(x=betweenness_value, y = letter, color = betweenness_type, size=2.5))  +
  geom_text(data = df_labels_betweenness, aes(x = midpoint, y = letter, label = paste0(round(diff_betweenness_percent,2), " %")), 
            vjust = 1.5, size = 3.5, color = "red") +
    geom_text(aes(label = betweenness_value, color = betweenness_type), vjust = -1, size=3.5) +  # Label the points with centrality_value
   scale_color_brewer(palette = "Set1", direction = 1) +
  scale_y_discrete(expand=c(0.05,0.05)) +
  theme_minimal() +
  labs(title = "Centrality analysis",
       subtitle= "Change of centrality based on sum edge betweenness in hospitals",
       y = " Hospitals",
       x = "Sum Edge Betweenness") +
  theme(plot.title=element_text(family ="bold", hjust=0.5),
        axis.text.y = element_text(size = 16),
        plot.subtitle = element_text(colour="#626262", hjust=0.5),
        legend.position='bottom') 

### table

df_hospitals_betweenness_both_compare <- hospitals_betweenness_both_long |> 
                                              mutate(diff_betweenness_percent = 
                                                          round((betwenness_post-betwenness_pre)/(betwenness_pre)*100,2),
                                                      diff = betwenness_post-betwenness_pre,
                                                      letter = LETTERS[1:22],
                                                      ds_cnes = stringr::str_to_title(ds_cnes))
                                                  
  
DT::datatable(subset(df_hospitals_betweenness_both_compare,
                     select=c("letter","cd_cnes","ds_cnes","betwenness_pre","betwenness_post","diff_betweenness_percent")),
              colnames=c("ID","Code","Name","Betweenness Pre","Betweenness Post","Diff(%)"),
              filter="top",
              extensions="Buttons",
                         options=list(
                           dom="Bfrtip",
                           buttons=c("copy","csv","pdf"),
                            initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'background-color': '#d50038', 'color': '#fff'});",
    "}")
                                  )
              ) |> 
    DT::formatStyle("betwenness_pre",
              background=DT::styleColorBar(range(df_hospitals_betweenness_both_compare$betwenness_pre),'#ee8b8b'),
                backgroundSize = '98% 88%',
  backgroundRepeat = 'no-repeat',
  backgroundPosition = 'center') |>
    DT::formatStyle("betwenness_post",
              background=DT::styleColorBar(range(df_hospitals_betweenness_both_compare$betwenness_pre),'#8ba1b3'),
                backgroundSize = '98% 88%',
  backgroundRepeat = 'no-repeat',
  backgroundPosition = 'center')

```

```{r}
#| eval: false
#| echo: false
#| code-summary: "R: Visualizing closeness results obtained with pgrouting"

## Import
hospitals_closeness_both <- st_read(eisenberg_connection, "hospitals_closeness_both") 
## pivot
hospitals_closeness_long <- hospitals_closeness_both |>
                                  mutate(diff_closeness_percent =
                                           (closeness_post-closeness_pre)/(closeness_pre)*100,
                                         diff = closeness_post-closeness_pre,
                                         letter = LETTERS[1:20]) |>
  arrange(letter) |> 
                             pivot_longer(
                                     c(closeness_pre,closeness_post),
                                      names_to="closeness_type",
                                      values_to="closeness_value") 
                   
hospitals_closeness_long_pre <- 
  hospitals_closeness_long |>
            filter(closeness_type == "closeness_pre")
hospitals_closeness_long_post <- 
  hospitals_closeness_long |>
  filter(closeness_type == "closeness_post")
##

df_labels <- hospitals_closeness_long |> 
  group_by(letter) |> 
  summarize(midpoint = mean(closeness_value), 
            diff_closeness_percent = first(diff_closeness_percent))
##
p <- ggplot(hospitals_closeness_long, aes(x=closeness_value, y=reorder(letter, desc(diff_closeness_percent)))) +
  geom_segment(data=hospitals_closeness_long_pre,
               aes(x=closeness_value, y = letter,
                   yend=hospitals_closeness_long_post$letter, xend=hospitals_closeness_long_post$closeness_value)) +
  geom_point(aes(x=closeness_value, y = letter, color = closeness_type, size=2.5))  +
  geom_text(data = df_labels, aes(x = midpoint, y = letter, label = paste0(round(diff_closeness_percent,2), " %")), 
            vjust = 1.5, size = 3.5, color = "red") +
    geom_text(aes(label = closeness_value, color = closeness_type), vjust = -1, size=3.5) +  # Label the points with centrality_value
   scale_color_brewer(palette = "Set1", direction = 1) +
  scale_y_discrete(expand=c(0.05,0.05)) +
  theme_minimal() +
  labs(title = "Centrality analysis",
       subtitle= "Change of centrality of hospitals due to the flooding",
       y = " Hospitals",
       x = "Closenness") +
  theme(plot.title=element_text(family ="bold", hjust=0.5),
        axis.text.y = element_text(size = 16),
        plot.subtitle = element_text(colour="#626262", hjust=0.5),
        legend.position='bottom') +
  xlim(0,47000) 

```

![](/media/accessibility_test2.png)

```{r}
#| code-summary: "R: Changes on the accesibiliy based on closeness"
hospitals_closeness_both <- st_read(eisenberg_connection, "hospitals_closeness_both") 
df_hospitals_closeness_both_compare <- hospitals_closeness_both |> 
                                              mutate(diff_closeness_percent = 
                                                          round((closeness_post-closeness_pre)/(closeness_pre)*100,2),
                                                      diff = closeness_post-closeness_pre,
                                                      letter = LETTERS[1:20],
                                                      ds_cnes = stringr::str_to_title(ds_cnes)) |> 
                                                  sf::st_drop_geometry() 
                                                  
  
DT::datatable(subset(df_hospitals_closeness_both_compare,
                     select=c("letter","cd_cnes","ds_cnes","closeness_pre","closeness_post","diff_closeness_percent")),
              colnames=c("ID","Code","Name","Closeness Pre","Closeness Post","Diff(%)"),
              filter="top",
              extensions="Buttons",
                         options=list(
                           dom="Bfrtip",
                           buttons=c("copy","csv","pdf"),
                            initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'background-color': '#d50038', 'color': '#fff'});",
    "}")
                                  )
              ) |> 
    DT::formatStyle("closeness_pre",
              background=DT::styleColorBar(range(df_hospitals_closeness_both_compare$closeness_pre),'#ee8b8b'),
                backgroundSize = '98% 88%',
  backgroundRepeat = 'no-repeat',
  backgroundPosition = 'center') |>
    DT::formatStyle("closeness_post",
              background=DT::styleColorBar(range(df_hospitals_closeness_both_compare$closeness_pre),'#8ba1b3'),
                backgroundSize = '98% 88%',
  backgroundRepeat = 'no-repeat',
  backgroundPosition = 'center')
```

### Closeness:

```{r}
#| code-summary: "R: Visualizing closeness of the healthcare facilities"

# Import he data
## hospital with closenesss values
closseness <-sf::st_read(eisenberg_connection, 
                            layer = "clossness_hospital_porto")
closseness_df <- closseness |>  arrange(ds_cnes, closeness) |>
                    mutate(lng= 
                              unlist(map(geom_hospital,1)),
                           lat=
                              unlist(map(geom_hospital,2)),
                           closeness_norm = 
                          (closseness$closeness - min(closseness$closeness)) / (max(closseness$closeness) - min(closseness$closeness)) * 100,
                          position = rank(-closeness))
## Create Color palette for visualization
pal <- colorQuantile(palette = "OrRd",closseness_df$closeness, n=4 )

## Create leaflet product

icons <- makeAwesomeIcon(
  icon = 'fa-heartbeat',
  iconColor = "#FFFFFF",
  markerColor = "#57142c",
  library = "fa"
)
leaflet(closseness_df) |>
    addProviderTiles(providers$OpenStreetMap.HOT) |>
    addCircles(data =closseness_df , radius = ~sqrt(closeness)*10, fillOpacity = .50, color =~pal(closeness)) |>
  addAwesomeMarkers(data=closseness_df,
                          icon =icons,
                          popup= ~paste0("<b>Código CNES: </b>", cd_cnes, "<br/>",
                                   "<b>Nome: </b>", ds_cnes, "<br/>",
                                   "<b> Closeness</b>:", closeness, "<br/>",
                                   "<b>Longitude: </b>", lng, "<br/>",
                                   "<b>Posição </b>", position, "<br/>"))

```

```{r}
#| code-summary: "R: Creating the table showing closeness values"
closseness <-sf::st_read(eisenberg_connection, 
                            layer = "clossness_hospital_porto")
closseness_no_geom <- closseness |>  
                          arrange(ds_cnes, closeness) |>
                          mutate(lng= 
                              unlist(map(geom_hospital,1)),
                           lat=
                              unlist(map(geom_hospital,2)),
                           closeness_norm = 
                          (closseness$closeness - min(closseness$closeness)) / (max(closseness$closeness) - min(closseness$closeness)) * 100,
                          position = rank(-closeness),
                          ds_cnes =stringr::str_to_title(ds_cnes)) |>
                  sf::st_drop_geometry()

DT::datatable(subset(closseness_no_geom, select=c("position","cd_cnes","ds_cnes","closeness")),
              extensions="Buttons",
                         options=list(
                           dom="Bfrtip",
                           buttons=c("copy","csv","pdf"),
                            initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'background-color': '#d50038', 'color': '#fff'});",
    "}")
                                  )
              ) |> 
    DT::formatStyle("closeness",
              background=DT::styleColorBar(range(closseness_no_geom$closeness),'#ee8b8b'),
                backgroundSize = '98% 88%',
  backgroundRepeat = 'no-repeat',
  backgroundPosition = 'center') 
```



## RQ3: Critical infrastructures

> Where are the most critical infrastructures located for accessing health facilities to reinforce urban resilience against flooding?