4  Extract difference

We extracted the difference between the road network and the flood extent creating the post-disaster network as result. In this section we described the process required to calculate the core metropolitan connectivity. For the intracity connectivity the process is the same but applied to each of the 9 networks.

We used a multi-step approach to extract the difference due to the demanding process in terms of computing resources. Similarly to the section detect flood, we also tested the naive approach extracting the difference at once creating spatial indexes in PostGIS or using duckdb (without spatial index¹.Given the results, we recommend the use of DuckDB for tasks with a high computational cost. However, due to the modest CPU used in this analysis, this multi-step was the feasible option.

The following map shows how we brokedown the process selecting several networks in each step. The layers pre-disaster network and post-disaster are the final results, while the networks roads in flood extent, boundary, outside and boundary-outside are intermediate steps.

# Select subset to quick visualization in the website
library(raster)
bounding_box_subset <-as(raster::extent(-51.248323,-51.092702,-30.048948,-29.947103),
                  "SpatialPolygons") 
proj4string(bounding_box_subset) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
bounding_box_subset <- bounding_box_subset |> sf::st_as_sfc()
# Flood extent 
flooding_cleaned_porto_union <- sf::st_read(connection,
                                    DBI::Id(
                                            schema="agile_gis_2025_rs",
                                            table="flooding_cleaned_porto_union"))
sf::sf_use_s2(FALSE) # This simplify the representation, we allowed it for visualization purpose.
flooding_cleaned_porto_union_subset <- sf::st_intersection(flooding_cleaned_porto_union,
                                                           bounding_box_subset)
## Read Road in flood extent
porto_ghs_road_in <- sf::st_read(connection,
                                    DBI::Id(
                                            schema="agile_gis_2025_rs",
                                            table="porto_ghs_road_in"))
## Subset
porto_ghs_road_in_subset <- sf::st_intersection(porto_ghs_road_in,
                                                bounding_box_subset) |> 
                                              rename(geometry=geom)
## Read Road out from flood extent
porto_ghs_road_out <- sf::st_read(connection,
                                    DBI::Id(
                                            schema="agile_gis_2025_rs",
                                            table="porto_ghs_road_out")) |>
                                            rename(geometry=geom)
## Subset
porto_ghs_road_out_subset <- sf::st_intersection(porto_ghs_road_out,
                                                 bounding_box_subset)
## Read Road on the flood extent boundary
porto_ghs_alegre_boundary <- sf::st_read(connection,
                                    DBI::Id(
                                            schema="agile_gis_2025_rs",
                                            table="porto_ghs_alegre_boundary"))
## Subset
porto_ghs_road_boundary_subset <- sf::st_intersection(porto_ghs_alegre_boundary,
                                                      bounding_box_subset) |>
                                                      rename(geometry=geom)
difference_outside_boundary <- sf::st_read(connection,
                                    DBI::Id(
                                            schema="agile_gis_2025_rs",
                                            table="difference_outside_boundary"))
## Subset
porto_ghs_road_outside_subset <- sf::st_intersection(difference_outside_boundary,
                                                     bounding_box_subset) |>
                                              rename(geometry=geom)
## Post & Pre disaster road network
porto_alegre_net_pre <- sf::st_read(connection,
                                    DBI::Id(
                                            schema="agile_gis_2025_rs",
                                            table="ghs_aoi_net_largest"))
## Subset
porto_ghs_road_net_pre_subset <- sf::st_intersection(porto_alegre_net_pre,
                                                     bounding_box_subset) |> 
                                                    rename(geometry = geom)
ghs_aoi_net_largest_post <-  sf::st_read(connection,
                                         DBI::Id(
                                                  schema="agile_gis_2025_rs",
                                                  table="ghs_aoi_net_largest_post")) |>
                                                      rename(geometry=geom)
## Subset
porto_ghs_road_net_post_subset <- sf::st_intersection(ghs_aoi_net_largest_post,
                                                      bounding_box_subset)

## Create dynamic map.
library(mapview)
library(leafpop)
mapview(porto_ghs_road_in_subset$geometry,
          color="red",
          lwd= 0.5,
          hide = TRUE,
          layer.name="A.Road in flood extent",
          popup=popupTable(porto_ghs_road_in_subset)) +
mapview(porto_ghs_road_outside_subset$geometry,
          color="lightgreen",
          lwd= 0.5,
          hide = TRUE,
          layer.name="D.Road boundary-outside",
          popup=popupTable(porto_ghs_road_outside_subset)) +
  mapview(porto_ghs_road_boundary_subset$geometry,
          color="yellow",
          lwd= 0.5,
          hide = TRUE,
          layer.name="C.Road boundary",
          popup=popupTable(porto_ghs_road_boundary_subset)) +    
  mapview(porto_ghs_road_out_subset,
          color="darkgreen",
          lwd= 0.2,
          hide = TRUE,
          layer.name="B.Road out from flood extent",
          popup=popupTable(porto_ghs_road_out_subset)) +     
mapview(porto_ghs_road_net_pre_subset,
          color="#a6611a",
          lwd= 0.2,
          layer.name="Pre-flooding network",
          popup=leafpop::popupTable(porto_ghs_road_net_pre_subset)) +
mapview(porto_ghs_road_net_post_subset,
          color="#018571",
          lwd= 0.2,
          layer.name="Post-flooding network",
          popup=leafpop::popupTable(porto_ghs_road_net_post_subset)) +    
  mapview(flooding_cleaned_porto_union_subset,
          color="darkblue",
          alpha.regions= 0.25,
          layer.name="Flood extent") +
  mapview(bounding_box_subset,
          color="red",
          alpha.regions = 0,
          lwd= 1,
          layer.name="Subset of the AoI") 

4.1 Multi-step approach

4.1.1 Road in flood extent

We selected the road segments in the flood extent applying the function st_intersection() with high computational cost only to the road segments inside the flood extent using the function st_contains(). The first two lines of the code [ORN-7] created spatial indexes to improve the performance.

CREATE INDEX idx_flooding_cleaned_porto_union_geom ON flooding_cleaned_porto_union USING gist(geom);
CREATE INDEX idx_ghs_aoi_net_largest ON ghs_aoi_net_largest USING gist(geom);
EXPLAIN ANALYZE
CREATE TABLE porto_ghs_road_in AS
  SELECT 
      net.id,
  CASE
    WHEN ST_Contains(flood.geom, net.geom)
    THEN net.geom
    ELSE st_intersection(net.geom, flood.geom)
  END AS geom
FROM
  ghs_aoi_net_largest AS net
JOIN
  flooding_cleaned_porto_union AS flood
ON st_intersects(net.geom, flood.geom);

---time:  11061.938 ms

4.1.2 Road out from flood extent

The code chunk [ORN-8] filtered out the road segments inside the flood by using a subquery of the table created in the code chunk [ORN-7]. Unlike the function st_intersection() in the previous code chunk that returned new geometries taking 11061 ms, using the numeric values in the where clause only took 359 ms.

--- 2) Network outside the flooding mask
EXPLAIN ANALYZE
CREATE TABLE porto_ghs_road_out AS
SELECT 
  net.*
FROM 
  ghs_aoi_net_largest AS net
WHERE 
  net.id NOT IN (SELECT 
                                id
                  FROM 
                        porto_ghs_road_in);
---time: 359.177 ms

4.1.3 Road on the flood extent boundary

We selected the road segments located in the boundary of the flood using the function st_exteriorring(). First, a CTE table named exterior_flood created a geometry that only contained the boundaries of the flood extent. Then a similar query to [ORN-7] selected the road network in the flood boundary.

CREATE INDEX idx_flooding_cleaned_porto ON flooding_cleaned_porto USING gist(geom);
EXPLAIN ANALYZE
CREATE TABLE porto_ghs_alegre_boundary AS
WITH exterior_flood AS (
SELECT
    ST_ExteriorRing((ST_dump(union_geom)).geom) AS geom
FROM
    (
    SELECT
        ST_Union(flood.geom) AS union_geom
    FROM
        flooding_cleaned_porto AS flood
    ) AS subquery
)
SELECT 
    net."id",
    CASE
        WHEN NOT ST_Contains(flood.geom, net.geom)
        THEN net.geom
        ELSE st_intersection(net.geom, flood.geom)
    END AS geom,
        net."toId",
        net."fromId",
        weight,
        "undrctI",
        "bdrctId"
FROM
    ghs_aoi_net_largest AS net
JOIN
    exterior_flood AS flood
    ON
    st_intersects(net.geom, flood.geom);
    
---time: 17600.814 ms

For the post-disaster network we were only interested in those segments that fell outside from the flood boundary. The code chunk [ORN-10] calculate the difference on the road segments found in the flood boundary. Although the function st_difference() has high computing costs, we only applied it to the small fraction of the road network located in the flood boundary.

CREATE INDEX idx_porto_ghs_alegre_boundary_geom ON porto_ghs_alegre_boundary USING gist(geom);
EXPLAIN ANALYZE
CREATE TABLE difference_outside_boundary AS
SELECT
    net.id,
    net."toId",
    net."fromId"",
    weight,
    net."undrctI",
    net."bdrctId",
    st_difference(net.geom, flood.geom) AS geom
FROM 
    porto_ghs_alegre_boundary AS net,
    flooding_cleaned_porto_union AS flood;
    
---time: 11235.65 ms  

The next query created the post-disaster network merging the road segments outside the flood extent with those located in the external side of the flood boundary. We included the name of the variables to ensure the same order, so the union is created sucessfuly.

CREATE INDEX idx_difference_outside_boundary ON difference_outside_boundary USING btree(id); 
CREATE INDEX idx_porto_ghs_road_out_id_ ON porto_ghs_road_out USING btree(id);
EXPLAIN ANALYZE
CREATE TABLE ghs_aoi_net_largest_post AS
SELECT
    id,
    geom,
    "toId",
    "fromId",
    weight,
    "undrctI",
    "bdrctId"
FROM
  porto_ghs_road_out
UNION
SELECT 
    id,
    geom,
    "toId",
    "fromId",
    weight,
    "undrctI",
    "bdrctId"
FROM 
  difference_outside_boundary;
---time: 581.093 ms

Troubleshooting 8 Important to use the united flood extent at the end, after using the subdivided for st_difference. Also mention the order of the columns Troubleshooting 9 When using the queries, a column is not found. Depending on how we imported the data, the name of the variables may change slightly. For example undrctI also is named unidirectid. Checking the names of the variables solved this problem

---

1. DuckDB: https://github.com/duckdb/duckdb-spatial