Show the code 
# Libraries
library(tidyverse)
library(terra)
library(tidyterra)
library(sf)
library(patchwork)
library(gt)
library(gtExtras)
library(purrr)

terraOptions(progress = 0)

3.1 Functions

In addition to the packages, we define our own functions for the analysis:

Show the code 
#' Function: get_perc_veg_buffer ----
#' Function for calculating if the percentage of vegetation within a buffer area is higher than a threshold 
#' @description It returns a raster indicating if the percentage of flammable vegetation within a circumference, from the centre of the grid (buffer), higher than a threshold (YES = 1, NO = 0)
#' @param .veg raster layer indicating if the vegetation is considered as fuel or not (YES = 1, NO = 0). 
#' @param .radio_buffer Radio of the circumference.
#' @param .threshold Proportion: from 0-1 (default  = 0.5; i.e., 50%) 
get_perc_veg_buffer <- function(.veg,
                                .radio_buffer,
                                .threshold = 0.5){
  
  # Percentage-vegetation 
  fwm <- terra::focalMat(.veg, d = radio_buffer, type = "circle", fillNA = TRUE)
  focal_perc <- terra::focal(.veg, w = fwm, fun = "sum", na.rm = TRUE) 
  
  # Higher than threshold (YES = 1, No = 0)
  wv_perc <- focal_perc
  wv_perc[wv_perc <  .threshold] <- 0
  wv_perc[wv_perc >= .threshold] <- 1
  
  return(wv_perc)
  
} 

#' Function: get_large_polygons ----
#' Function for removing small areas from a raster file. Used for calculating large forest areas (e.g., > 5 km2).  
#' @description It returns a raster with all the polygons larger than a minimal area (YES = 1, NO = 0)
#' @param .x starts object (i.e., 1 - fuel vegetation, 0 - no fuel vegetation). 
#' @param .min_area Minimum area [in km2] for considering an vegetated area as large  (Default 5 km^2)
#' @param .resolution Resolution of the final raster (default - 100 x 100 m)
get_large_polygons <- function(.x,               
                               .min_area = 5,    
                               .resolution = 100) {
  
  df_v <- .x |> 
    # Convert to polygons
    st_as_sf(as_points = FALSE, merge = TRUE) %>% 
    # Calculate the area of the polygons
    mutate(area = units::set_units(st_area(.), "km2")) |> 
    # Remove lower than 5 km2
    filter(area > units::set_units(.min_area, "km2")) |> 
    # Set columns names
    setNames(c("large_wild_veg", "geometry", "area")) |> 
    # Remove area column
    select(-area) |> 
    # transform to SpatVector  
    vect() |> 
    # aggregate in one polygon
    aggregate(by = "large_wild_veg") |> 
    select(large_wild_veg)
  
  # Rasterize
  df_geom <- rast(df_v, resolution = .resolution)
  df_r <- rasterize(df_v, df_geom)
  
  return(df_r)
  
}

#' get_dist_max ----
#' Get cells within a distance to vegetated area (YES --> 1; NO --> 0)  
#' @description It returns a raster with indicating the cells that are within the distance to the vegetated areas (YES = 1, NO = 0)
#' @param .x Raster with distances (from `terra::distance`) 
#' @param .max_distance Maximum distance in [m] (e.g., 2400 m)
get_dist_max <- function(.x, .max_dist) {
  
  dist <- .x
  dist[dist <= .max_dist] <- 1
  dist[dist >  .max_dist] <- 0
  
  return(dist)
  
}

#' Function: get_wui_raster ----
#' Function for WUI classification
#' @description it returns a raster with WUI areas (Interface or Intermix)
#' @param .veg_perc raster layer indicating if the percentage of vegetation in
#' the buffer area is higher than a threshold (e.g., > 50%). (YES = 1, NO = 0). 
#' @param .veg_dist raster layer indicating if the a grid cell is within a 
#' distance (e.g., 2.4 km) of a large area of flammable vegetation (YES = 1, NO = 0)
#' @param .pot_wui raster layer with potential WUI grids based on building density (`wui_potential`).
#' @param .region SpatVector layer with the region area (`region_v`) 
get_wui_raster <- function(.veg_perc,
                           .veg_dist,
                           .pot_wui = wui_potential,
                           .region = region_v){

  # Potential WUI grid ----
  wui_potential <- .pot_wui
  wui_potential <- as.factor(wui_potential)
  
  # Match extents
  perc <- resample(.veg_perc, wui_potential, method = "mode")
  dist <- resample(.veg_dist, wui_potential, method = "mode")
  
  # Intermix 
  wui <- perc + wui_potential # 2 = WUI intermix
  wui[wui == 1] <- 0          # 1 = Potential interface --> 0
  
  # Interface 
  wui <- wui + dist
  wui[wui == 0] <- NA # No WUI
  wui[wui == 1] <- 1  # Interface
  wui[wui == 2] <- 2  # Intermix
  wui[wui == 3] <- 2  # Intermix
  
  # NAs as 0 inside the study area (otherwise NA)
  wui <- ifel(is.na(wui), 0, wui)
  
  # Mask to region (e.g., NA = ourside norway)
  wui <- terra::mask(wui, .region)
  
  # Identify levels (1 = Interface, 2 = Intermix)
  levels(wui) <- data.frame(id = 0:2,
                            wui_class = c("No_wui", "Interface", "Intermix")
                            )
  
  # Add colour table 
  coltab(wui) <- data.frame(id = 0:2, 
                            col = c("#FAFAFA", "#FF6633", "#003366"))
  
  
  # Return SpatRaster
  set.values(wui)
  return(wui)
  
}

#' add_colours_wui 
#' Function for adding WUI colours in a SpatRaster 
#' @description Add colours to a WUI raster layer (No_WUI, Interface, Intermix)
#' @param .veg_perc raster layer indicating if the percentage of vegetation in
add_colours_wui <- function(.x) {
  
  r <- .x
  # Identify levels (1 = Interface, 2 = Intermix)
  df_lev <- data.frame(id = 0:2, wui_class = c("No_wui", "Interface", "Intermix"))
  df_lev_lst <- replicate(nlyr(r), df_lev, simplify = FALSE)
  levels(r) <- df_lev_lst
  # Add colour table 
  df_col <- data.frame(id = 0:2, col = c("#FAFAFA", "#FF6633", "#003366"))
  df_col_lst <- replicate(nlyr(r), df_col, simplify = FALSE)
  coltab(r) <- df_col_lst
  
  return(r)
  
}

3.2 Data preprocessing

3.2.1 Parameter

Show the code 
raster_resolution <- 250 # m
radio_buffer <- 500 # m
area_buffer <- pi * (radio_buffer/1000)^2 # km2

At national scale, we have performed a raster analysis by grids cells of 250 m x 250 m.

3.2.2 Study area

We load the last release of the county data from GISCO (year = 2021), with a resolution of 1:1 million and crs= 3035 (Figure 3.1). Due to size of the file, we perform the analysis by tiles.

Show the code 
# Create local directory for caching big files
dir.create("C:/GISCO_spatial_data")
options(gisco_cache_dir = "C:/GISCO_spatial_data")

# Counties ----
country_list <- c("NO", "SE")
counties <- giscoR::gisco_get_nuts(country = country_list,
                                   year = "2021",
                                   nuts_level = 3,
                                   epsg = "3035",
                                   resolution = "01") |> 
  filter(NUTS_NAME != "Svalbard") |> 
  filter(NUTS_NAME != "Jan Mayen") |> 
  select(NUTS_ID, NUTS_NAME, CNTR_CODE)

# Study region (SF object) ----
region_sf <- counties |> 
  summarize()

country_sf <- counties |>
  group_by(CNTR_CODE) |> 
  summarize()

# # Rogaland ----
# rogaland <- counties |> 
#   filter(NUTS_NAME == "Rogaland")

# Transform to terra::SpatVector (region_v) and SpatRaster (region_r)) -----
region_v <- vect(region_sf)
geometry_raster <- rast(region_v, resolution = raster_resolution)
region_r <- rasterize(region_v, geometry_raster)

no_v <- country_sf |> 
  filter(CNTR_CODE == "NO") |> 
  vect()

se_v <- country_sf |> 
  filter(CNTR_CODE == "SE") |> 
  vect()

# Tiles 
raster_tiles <- rast(nrows = 15,
                     ncols = 15,
                     ext(region_r),
                     crs = "EPSG:3035")

dir.create("data/processed_data")
dir.create("data/processed_data/tiles")
dir.create("data/processed_data/tiles/region")

region_tiles <- makeTiles(
  x = region_r,
  y = raster_tiles,
  filename = "data/processed_data/tiles/region/region_r_.tif",
  na.rm = TRUE,
  extend = TRUE
)

region_vrt <- vrt(region_tiles)

# Plot Study area ----

# world 
world <- rworldmap::getMap(resolution = "high") |> 
  st_as_sf() |> 
  st_transform(crs = 3035) 

p1 <- ggplot() +
  geom_sf(data = world) +
  geom_sf(data = country_sf, fill = c("#D55E00", "#4E84C4")) +
  coord_sf(xlim = c(2500000, 5402285),
           ylim = c(1250000, 5500000),
           expand = FALSE) +
  theme_bw() +
  theme(panel.background = element_rect(fill = "lightgrey")) +
  theme(panel.grid = element_blank(),
        axis.text  = element_blank(),
        axis.ticks = element_blank()) 
  
p2 <- ggplot() +
  geom_sf(data = region_sf, fill = "#52854C") +
  geom_sf(data = country_sf, fill = NA) +
  theme_bw() 

p2 + inset_element(p1, 0.03, 0.70, 0.6, 0.95, align_to = 'full')
Figure 3.1: Study area (Norway and Sweden)

3.2.3 Building density

We have downloaded the building data from OpenStreetMaps (geofabrik: link), in pbf format. The version of the datasets are:

Name Last modified Size
norway-latest-osm.pbf 2023-08-17 04:37 1.2G
sweden-latest-osm.pbf 2023-11-16 04:50 656M

They are big files and we have not put them on GitHub. If you wanted to replicate the document, you would need to download the file and save it with the same path as here, or change the path in the code (NOTE: you may have small differences in the results if OSM updates the data). To calculate the building density we do the following steps:

  • Transform the .pbf file to .gpkg, extracting only multipolygons (run only one time)
  • Read building polygons
  • Calculate the centroid of each building, assigning a value of 1 to the building (building = 1)
  • Calculate the building density in a grid cell of 250 m x 250 m (it is done by tiles). It is the sum of all the buildings (building = 1) within a grid cell
  • Merge all tiles in one raster (terra::mosaic)
  • Smooth the data with a circular moving window of 500 m; i.e., buffer area (terra::focal)
  • The density building of a grid cell is therefore the number of buildings in the buffer area divided by the area of the buffer. Potential WUI areas are those grid cells with a building density higher than 6.17 #/km^2
Show the code 
# # Transform OSM data to gpkg (only buildings) -- RUN ONLY ONCE 
# osm_no_pbf <- "data/big_data/OSM/norway-latest.osm.pbf"
# osm_no_building_gpkg <- osmextract::oe_vectortranslate(osm_pbf,
#                                                 layer = "multipolygons")

# # Transform OSM data to gpkg (only buildings) -- RUN ONLY ONCE
# osm_se_pbf <- "data/big_data/OSM/sweden-latest.osm.pbf"
# osm_se_building_gpkg <- osmextract::oe_vectortranslate(osm_se_pbf,
#                                                 layer = "multipolygons")

#' Function for getting buildings counts from OSM
#' @description Read only buildings from OSM.gpkg file (by tiles) and calculate building density
#' @param .tile_path Tiles where building counts should be calculated 
#' @param .data_path path to folder with OSM data (.gpkg format) 
get_building_counts <- function(.tile_path,
                                .data_path = "data/big_data/OSM"){ 
  
  # read tile
  tile <- terra::rast(.tile_path)
  files <- list.files(.data_path,
                      pattern="*.gpkg",
                      full.names = TRUE)
  
  # Get bounding box tile
  wkt_tile <- tile |>
    terra::project("EPSG:4326") |> 
    terra::ext() |> 
    terra::vect() |>
    terra::geom(wkt = TRUE) 
  
  # Calculate number buildings
  bd <- lapply(files, function(.x) { sf::read_sf(.x, wkt_filter = wkt_tile) }) |>  
    dplyr::bind_rows() |> 
    tidyr::drop_na(building) |> 
    dplyr::select(building) |>
    sf::st_transform(crs = 3035) |> 
    sf::st_centroid() |> 
    dplyr::mutate(building = 1) |> 
    terra::vect() |> 
    terra::rasterize(tile, fun = sum) 
  
  return(bd)
  
}

# Building density ----
build_counts <- region_tiles %>%
  # Density by tiles
  map(get_building_counts) |> 
  # Merge tiles
  sprc() |> 
  mosaic(fun = "mean")
terra::writeRaster(build_counts, 
                   "data/processed_data/build_counts.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Focal (Smooth values) ----
# Number of buildings in a circle
fwm <- focalMat(build_counts, d = radio_buffer, type = "circle", fillNA = TRUE)
fwm[fwm > 0] <- 1
build_focal <- focal(build_counts, w = fwm, fun = sum , na.rm = TRUE)
terra::writeRaster(build_focal, 
                   "data/processed_data/build_focal.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)
Show the code 
# Read data
build_counts <- terra::rast("data/processed_data/build_counts.tif")
build_focal <- terra::rast("data/processed_data/build_focal.tif")

# Plot ----
p1 <- ggplot() +
  tidyterra::geom_spatvector(data = region_v, fill = "white") +
  tidyterra::geom_spatraster(data = build_counts) +
  scale_fill_whitebox_c(name = "N. Buildings",
                        palette = "viridi",
                        trans = "log10") +
  labs(title = "Original data from OSM") +
  theme_void()

p2 <- ggplot() +
  tidyterra::geom_spatvector(data = region_v, fill = "white") +
  tidyterra::geom_spatraster(data = build_focal) +
  scale_fill_whitebox_c(name = "N. Buildings",
                        palette = "viridi",
                        trans = "log10") +
   labs(title = "Focal values (circular window)") +
  theme_void()

p1 + p2
Figure 3.2: Number of buildings per grid cells of 250 m x 250 m

3.2.4 Fuel/vegetation maps

3.2.4.1 EFFIS Fuel Map

The EFFIS Fuel Map (EFFIS (2017)) is a pan-European map that classifies the vegetation in 10 categories according to the fire behaviour models developed by Anderson, 1982 (Anderson (1982)) (i.e., it excludes the last 3 categorise of the model). Resolution 250 x 250 m.

We remove small polygons (lower than 5 km2), and calculate the distance from each grid cell to the nearest large wild vegetation. Then, we calculate the grids that are within 2.4 km (value = 1) or not (value = 0).

We also calculate the percentage of vegetation in the buffer area (*_wv_per). In the study region, there are only 8 EFFIS categories (Figure 3.3).

Show the code 
# # Dowload and unzip data (RUN once)
# effis_fuel_url <- "https://effis-gwis-cms.s3-eu-west-1.amazonaws.com/effis/applications/data-and-services/FuelMap_LAEA.zip"
# dir.create("data/big_data/EFFIS")
# download.file(effis_fuel_url, "data/big_data/EFFIS/FuelMap.zip")
# unzip(zipfile = "data/big_data/EFFIS/FuelMap.zip",
#       exdir = "data/big_data/EFFIS")

# Forest layer norway
effis_fuel_path <- "data/big_data/EFFIS/FuelMap2000_NFFL_LAEA.tif"
effis_fuel <- rast(effis_fuel_path) |> 
  crop(region_v) |> 
  mask(region_v)

# Wildland vegetation
effis_wild_veg <- effis_fuel
effis_wild_veg[effis_wild_veg >= 1] <- 1
names(effis_wild_veg) <- "wild_veg"
terra::writeRaster(effis_wild_veg, 
                   "data/processed_data/effis_wild_veg.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE) 

# Large forest areas (> 5km2)
effis_wild_veg_large <- stars::read_stars("data/processed_data/effis_wild_veg.tif") |> 
  get_large_polygons() 
terra::writeRaster(effis_wild_veg_large, 
                   "data/processed_data/effis_wild_veg_large.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Distances to large vegetated areas ----
effis_wild_veg_large_dist <- terra::distance(effis_wild_veg_large)
terra::writeRaster(effis_wild_veg_large_dist, 
                   "data/processed_data/effis_wild_veg_large_dist.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Percentage-vegetation 
effis_wv_perc <- get_perc_veg_buffer(
  .veg = effis_wild_veg,
  .radio_buffer = radio_buffer)
terra::writeRaster(effis_wv_perc, 
                   "data/processed_data/effis_wv_perc.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Plot ----
effis_wild_veg_list_no <- tribble(~code, ~cover_class, ~RGB,
                               1, "Short grass",    "16-240-96",
                               3, "Tall grass", "16-128-16",
                               5, "Brush",  "255-255-0",
                               6, "Dormant brush", "255-204-0",
                               7, "Southern rough", "255-102-0",
                               8, "Closed timber litter", "192-128-64",
                               9, "Hardwood litter", "160-64-0",
                               10, "Timber", "112-48-48") |>
  separate(RGB, c("red", "green", "blues"), sep = "-") |>  
  mutate(rgb_code = rgb(red, green, blues, maxColorValue = 255))

# Plot
effis_fuel_cls <- effis_fuel
cls <- data.frame(id = effis_wild_veg_list_no$code, 
                  value = effis_wild_veg_list_no$cover_class)
levels(effis_fuel_cls) <- cls

plot(region_r, col = "lightgrey", axes = FALSE, legend = FALSE)
plot(effis_fuel_cls, 
     col = effis_wild_veg_list_no$rgb_code,
     axes = FALSE, add = T)
Figure 3.3: Wildland vegetation in Norway by type of vegetation (Source EFFIS Fuel Map)
Show the code 
# Read Effis data
effis_wild_veg <- terra::rast("data/processed_data/effis_wild_veg.tif")

# Get cells within 2400 m (YES --> 1; NO --> 0)
dist_path <- "data/processed_data/effis_wild_veg_large_dist.tif"
effis_wild_veg_large_dist_2400 <- terra::rast(dist_path) |> 
  get_dist_max(2400)

# Plot
par(mfrow = c(1,2), adj = 0)

# Wildand vegetation (Yes = 1, No = NA)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(effis_wild_veg, col = "#009E73", axes = F, legend = F, add = T)
title("A)", adj = 0)
# Within a distance of 2.4 km (Yes = 1, No = 0)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(effis_wild_veg_large_dist_2400, col = c(NA, "#0072B2"), axes = F, legend = F, add = T)
title("B)", adj = 0)
Figure 3.4: a) Wildland vegetation in Norway, and b) grid cells within a distance to large vegetation areas (>5km2) of 2.4 km (Source EFFIS Fuel Map)

3.2.4.2 Corine Land Cover

We have used the Corine Land Cover 2018 (Version 2020_20u1), downloaded from Copernicus. The resolution of the raster is 100 m x 100 m. We classify as potential wildland vegetation the following classes Figure 3.5:

Show the code 
# Load data ----
clc2018_path <- "data/big_data/corine/DATA/U2018_CLC2018_V2020_20u1.tif" 

clc2018_r <- rast(clc2018_path) |> 
  crop(region_v) |> 
  mask(region_v)

# Get Wildland vegetation ----
clc2018_wild_veg_list <- tribble(~CODE_18, ~LABEL3, ~RGB,
                                 311,   "Broad-leaved forest",  "128-255-000",
                                 312,   "Coniferous forest",    "000-166-000",
                                 313,   "Mixed forest", "077-255-000",
                                 321,   "Natural grasslands",   "204-242-077",
                                 322,   "Moors and heathland",  "166-255-128",
                                 323,   "Sclerophyllous vegetation", "166-230-077",
                                 324,   "Transitional woodland-shrub",  "166-242-000") |>
  separate(RGB, c("red", "green", "blues"), sep = "-") |>  
  mutate(rgb_code = rgb(red, green, blues, maxColorValue = 255))

# Subset target vegetation
corine_fuel_veg <- clc2018_r |> 
  terra::subst(clc2018_wild_veg_list$LABEL3, clc2018_wild_veg_list$LABEL3, others = NA) |> 
  droplevels() 

# Wildland vegetation: yes --> 1 
corine_wild_veg <- as.int(corine_fuel_veg)
corine_wild_veg[corine_wild_veg > 0] <- 1
# Remove colours associated with the rasters
coltab(corine_wild_veg) <- NULL
names(corine_wild_veg) <- "wild_veg"
terra::writeRaster(corine_wild_veg, 
                   "data/processed_data/corine_wild_veg.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Large forest areas (> 5km2)
corine_wild_veg_large <- stars::read_stars("data/processed_data/corine_wild_veg.tif") |> 
  get_large_polygons() 
terra::writeRaster(corine_wild_veg_large, 
                   "data/processed_data/corine_wild_veg_large.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Distances to large vegetated areas ----
set.values(corine_wild_veg_large)
corine_wild_veg_large_dist <- terra::distance(corine_wild_veg_large)
terra::writeRaster(corine_wild_veg_large_dist, 
                   "data/processed_data/corine_wild_veg_large_dist.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Percentage-vegetation ----
corine_wv_perc <- get_perc_veg_buffer(
  .veg = corine_wild_veg,
  .radio_buffer = radio_buffer)
terra::writeRaster(corine_wv_perc, 
                   "data/processed_data/corine_wv_perc.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Plot ----
plot(region_r, col = "lightgrey", axes = FALSE, legend = FALSE)
plot(corine_fuel_veg, col = clc2018_wild_veg_list$rgb_code, 
     axes = FALSE, add = T)
Figure 3.5: Wildland vegetation in Norway by type of vegetation (Source CORINE)
Show the code 
# Read CORINE data
corine_wild_veg <- terra::rast("data/processed_data/corine_wild_veg.tif")

# Get cells within 2400 m (YES --> 1; NO --> 0)
dist_path <- "data/processed_data/corine_wild_veg_large_dist.tif"
corine_wild_veg_large_dist_2400 <- terra::rast(dist_path) |> 
  get_dist_max(2400)


# Plot
par(mfrow = c(1,2), adj = 0)

# Wildland vegetation (Yes = 1, No = NA)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(corine_wild_veg, col = "#009E73", axes = F, legend = F, add = T)
title("A)", adj = 0)
# Within a distance of 2.4 km (Yes = 1, No = 0)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(corine_wild_veg_large_dist_2400, col = c(NA, "#0072B2"), axes = F, legend = F, add = T)
title("B)", adj = 0)
Figure 3.6: a) Wildland vegetation in Norway, and b) grid cells within a distance to large vegetation areas (>5km2) of 2.4 km (Source CORINE)

3.2.4.3 National datasets

3.2.4.3.1 Norwegian Forest resource map - SR16

We downloaded the Norwegian Forest Resource map (Skogressurskart: SR16 - raster format: resolution 16 m) from NIBIO. We downloaded the raster data by county (Fylke) in EUREF89 UTM 33, 2d (25833). There are multiple raster files but we only used the dominant tree species (SRTRESLAG), which is classified in three categories:

- gran - spruce: 1

- furu - pine: 2

- lauv - deciduous: 3

We consider all of them as wildland vegetation (we reclassified the raster as 1 = wildland vegetation, and 0 = No wildland vegetation).

We aggregate the data by a factor of 6.25, so we past from an original resolution of 16x16 meters to 100x100 meters.

Show the code 
# Vegetation layer norway ----
sr16r_path <- list.files("data/big_data/sr16/srtreslag",
                         pattern = "*.tif",
                         full.names = TRUE)
sr16_r <- map(sr16r_path, rast)

# Wildland vegetation
get_wild_veg <- function(.x){
    r <- .x
    r[r >= 1] <- 1 
    return(r)
  }
sr16_wild_veg <-  map(sr16_r, get_wild_veg) 
  
# Aggregate 
change_resolution <- function(.x, .fact = 6.25) {
  aggregate(.x, fact = .fact, fun = "modal") 
}

sr16_wild_veg <- map(sr16_wild_veg, change_resolution)

# Merge rasters
sr16_wild_veg <- sr16_wild_veg |> 
  sprc() |> 
  mosaic(fun = "max") |> 
  project("EPSG:3035", method = "near")

# Export as .tif
set.values(sr16_wild_veg)
names(sr16_wild_veg) <- "wild_veg"
terra::writeRaster(sr16_wild_veg,
                   "data/processed_data/sr16_wild_veg.tif", 
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Large forest areas (> 5km2)
sr16_wild_veg_large <- stars::read_stars("data/processed_data/sr16_wild_veg.tif") |>   get_large_polygons() 
set.values(sr16_wild_veg_large)
terra::writeRaster(sr16_wild_veg_large, 
                   "data/processed_data/sr16_wild_veg_large.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Distances to large vegetated areas ----
sr16_wild_veg_large_dist <- terra::distance(sr16_wild_veg_large)
set.values(sr16_wild_veg_large_dist)
terra::writeRaster(sr16_wild_veg_large_dist, 
                   "data/processed_data/sr16_wild_veg_large_dist.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Percentage-vegetation ----
sr16_wv_perc <- get_perc_veg_buffer(
  .veg = sr16_wild_veg,
  .radio_buffer = radio_buffer)
terra::writeRaster(sr16_wv_perc, 
                   "data/processed_data/sr16_wv_perc.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)
Show the code 
# Plot original data ----
categorical_levels <- function(.x){
  r <- .x
  levels(r) <- data.frame(id = 1:3, type =  c("spruce", "pine", "deciduous"))
  return(r)
}

sr16_vrt <- terra::vrt(sr16r_path) 
levels(sr16_vrt) <- data.frame(id = 1:3, type =  c("spruce", "pine", "deciduous"))
no_v_p <- project(no_v, "EPSG:25833")
plot(no_v_p, col = "lightgrey", axes = FALSE, legend = FALSE)
plot(sr16_vrt, col = c("lightgreen", "darkgreen", "yellow"), axes = FALSE, add = T)
Figure 3.7: Wildland vegetation in Norway by type of vegetation (Source SR16)
3.2.4.3.2 Swedish National Land Cover Database (NMD)

We download the NMD Base map from the Swedish Environmental Protection Agency (Nationella marktäckedata 2018; basskikt). It classified Sweden in 6 main levels, which later are subdivided in level 2 and level 3 (Figure 3.8):

NMD categories (table Bilaga 1: Nomenklatur)
Level 1 Level 2 Level 3
  1. Skog
 1.1 Skog utanför våtmark

1.1.1 Tallskog utanför våtmark

1.1.2 Granskog utanför våtmark

1.1.3 Barrblandskog utanför våtmark

1.1.4 Lövblandad barrskog utanför våtmark

1.1.5 Triviallövskog utanför våtmark

1.1.6 Ädellövskog utanför våtmark

1.1.7 Triviallövskog med ädellövinslag utanför våtmark

1.1.8 Temporärt ej skog utanför våtmark
1.2 Skog på våtmark

1.2.1 Tallskog på våtmark

1.2.2 Granskog på våtmark

1.2.3 Barrblandskog på våtmark

1.2.4 Lövblandad barrskog på våtmark

1.2.5 Triviallövskog på våtmark

1.2.6 Ädellövskog på våtmark

1.2.7 Triviallövskog med ädellövinslag på våtmark

1.2.8 Temporärt ej skog på våtmark
  1. Våtmark
  1. Åkermark
  1. Övrig öppen mark
 4.1 Övrig öppen mark utan vegetation
4.2 Övrig öppen mark med vegetation
  1. Exploaterad mark
 5.1 Exploaterad mark, byggnad
5.2 Exploaterad mark, ej byggnad eller väg/järnväg
5.3 Exploaterad mark, väg/järnväg
  1. Vatten
 6.1 Sjö och vattendrag
6.2 Hav
Show the code 
nmd_path <- "data/big_data/NMD2018/nmd2018bas_ogeneraliserad_v1_1.tif" 
nmd_r <- rast(nmd_path) 

# Vegetation code
levles_lst <- levels(nmd_r)
levles_lst[[1]]$Klass <- iconv(levles_lst[[1]]$Klass, "latin1")
levels(nmd_r) <- levles_lst

nmd_veg_list <- levles_lst[[1]] |>
  as_tibble() |> 
  mutate(Klass = ifelse(Klass == "", NA, Klass)) |> 
  drop_na()

plot(nmd_r)
Figure 3.8: Vegetation map of Sweden (National Land Cover Database 2018)

Then, we classified as non-burnable the following categories:

Show the code 
nmd_no_fuel_list <- nmd_veg_list |> 
  filter(value %in% c(3, 41, 42, 51, 52, 53, 61, 62, 131, 169, 170, 179, 181, 189, 190))

nmd_no_fuel_list |> 
  gt()
Table 3.1: Non-burnable vegetation (categories)
value Klass
3 Åkermark
41 Övrig öppen mark utan vegetation
42 Övrig öppen mark med vegetation
51 Exploaterad mark, byggnad
52 Exploaterad mark, ej byggnad eller väg/järnväg
53 Exploaterad mark, väg/järnväg
61 Sjö och vattendrag
62 Hav
131 Åkermark
169 Övrig öppen mark utan vegetation
170 Övrig öppen mark med vegetation
179 Exploaterad mark, byggnad
181 Exploaterad mark, väg/järnväg
189 Sjö och vattendrag
190 Hav

The original resolution of the map is 10 meters. We aggregated it (by a factor of 10) and transformed the coordinate system from SWEREF99 TM (EPSG:3006) to ETRS89-extended / LAEA Europe (EPSG:3035).

Show the code 
# Reduce resolution 
nmd_r_a <- aggregate(nmd_r, fact = 10, fun = "modal") |> 
  project("EPSG:3035", method = "near")

# Remove no fuel vegetations 
nmd_wild_veg <- nmd_r_a |> 
  mutate(Klass = ifelse(Klass %in% nmd_no_fuel_list$Klass, NA, Klass))

# Wild vegetation; YES = 1, NO = 0
nmd_wild_veg[nmd_wild_veg == 1] <- NA  # Outside Sweden 
nmd_wild_veg[nmd_wild_veg  > 1] <- 1

# Export as .tif
names(nmd_wild_veg) <- "wild_veg"
coltab(nmd_wild_veg) <- NULL
set.values(nmd_wild_veg)
terra::writeRaster(nmd_wild_veg,
                   "data/processed_data/nmd_wild_veg.tif", 
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Large forest areas (> 5km2)
nmd_wild_veg_large <- stars::read_stars("data/processed_data/nmd_wild_veg.tif") |> 
  get_large_polygons() 
set.values(nmd_wild_veg_large)
terra::writeRaster(nmd_wild_veg_large, 
                   "data/processed_data/nmd_wild_veg_large.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Distances to large vegetated areas ----
nmd_wild_veg_large_dist <- terra::distance(nmd_wild_veg_large)
set.values(nmd_wild_veg_large_dist)
terra::writeRaster(nmd_wild_veg_large_dist, 
                   "data/processed_data/nmd_wild_veg_large_dist.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Percentage-vegetation ----
nmd_wv_perc <- get_perc_veg_buffer(
  .veg = nmd_wild_veg,
  .radio_buffer = radio_buffer)
set.values(nmd_wv_perc)
terra::writeRaster(nmd_wv_perc, 
                   "data/processed_data/nmd_wv_perc.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

Finally, we combined both maps in one (Figure 3.9)

Show the code 
# Read SR16 data
sr16_wild_veg <- terra::rast("data/processed_data/sr16_wild_veg.tif")
dist_path <- "data/processed_data/sr16_wild_veg_large_dist.tif"
sr16_wild_veg_large_dist_2400 <- terra::rast(dist_path) |> 
  get_dist_max(2400)

# Read NRM data
nmd_wild_veg <- terra::rast("data/processed_data/nmd_wild_veg.tif")
dist_path <- "data/processed_data/nmd_wild_veg_large_dist.tif"
nmd_wild_veg_large_dist_2400 <- terra::rast(dist_path) |> 
  get_dist_max(2400)

# Plot
par(mfrow = c(1,2), adj = 0)

# Wildland vegetation (Yes = 1, No = NA)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(sr16_wild_veg, col = "#009E73", axes = F, legend = F, add = T)
plot(nmd_wild_veg, col = "#009E73", axes = F, legend = F, add = T)
title("A)", adj = 0)
# Within a distance of 2.4 km (Yes = 1, No = 0)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(sr16_wild_veg_large_dist_2400, 
     col = c(NA, "#0072B2"),
     axes = F, 
     legend = F,
     add = T)
plot(nmd_wild_veg_large_dist_2400, 
     col = c(NA, "#0072B2"), 
     axes = F, 
     legend = F, 
     add = T)
title("B)", adj = 0)
Figure 3.9: a) Wildland vegetation in the region, and b) grid cells within a distance to large vegetation areas (>5km2) of 2.4 km (Source SR16 and NMD)

3.2.4.4 FirEUrisk_Europe_fuel_map

We use the fuel map produced by (Aragoneses, Garcia, and Chuvieco 2022). The original resolution of the map is 1km, so we resampled it to a grid cell of 250x250 metres (same grids than region_r) based in the type that is most presented in the area; i.e., mode (Figure 3.10).

Show the code 
# Codes
firEUrisk_fuel_type <- tribble(~code, ~description,
                                1111,"Open broadleaf evergreen forest",
                                  23, "High shrubland",
                                1112, "Closed broadleaf evergreen forest",
                                  31, "Low grassland",
                                1121, "Open broadleaf deciduous forest",
                                  32, "Medium grassland",
                                1122, "Closed broadleaf deciduous forest",
                                  33, "High grassland",
                                1211, "Open needleleaf evergreen forest",
                                  41, "Herbaceous cropland",
                                1212, "Closed needleleaf evergreen forest",
                                  42, "Woody cropland",
                                1221, "Open needleleaf deciduous forest",
                                  51, "Wet and peat/semi-peat land – tree",
                                1222, "Closed needleleaf deciduous forest", 
                                  52, "Wet and peat/semi-peat land – shrubland",
                                1301, "Open mixed forest",
                                  53, "Wet and peat/semi-peat land – grassland",
                                1302, "Closed mixed forest",
                                  61, "Urban continuous fabric",
                                  21, "Low shrubland",
                                  62, "Urban discontinuous fabric",
                                  22, "Medium shrubland",
                                   7, "Nonfuel")

firEUrisk_no_fuel_type <- tribble(~code, ~description,
                                     41, "Herbaceous cropland",
                                     42, "Woody cropland",
                                     51, "Wet and peat/semi-peat land – tree",
                                     52, "Wet and peat/semi-peat land – shrubland",
                                     53, "Wet and peat/semi-peat land – grassland",
                                     61, "Urban continuous fabric",
                                     62, "Urban discontinuous fabric",
                                      7, "Nonfuel")

# Forest layer norway
firEUrisk_fuel_path <- "data/big_data/firEUrisk/FirEUrisk_Europe_fuel_map.tif"
firEUrisk_fuel <- rast(firEUrisk_fuel_path) |> 
  crop(region_v) |> 
  mask(region_v)

# Resample to grid cells of 250m x 250m (same than region_r)
firEUrisk_fuel_250m <- resample(firEUrisk_fuel, region_r, method = "mode") |> 
 # Remove no fuel vegetations 
  mutate(FirEUrisk_Europe_fuel_map = ifelse(FirEUrisk_Europe_fuel_map %in% firEUrisk_no_fuel_type$code,
                    NA,
                    FirEUrisk_Europe_fuel_map))

# Wildland vegetation
firEUrisk_wild_veg <- firEUrisk_fuel_250m 
firEUrisk_wild_veg[firEUrisk_wild_veg >= 1] <- 1
names(firEUrisk_wild_veg) <- "wild_veg"
terra::writeRaster(firEUrisk_wild_veg, 
                   "data/processed_data/firEUrisk_wild_veg.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Large forest areas (> 5km2)
firEUrisk_wild_veg_large <- stars::read_stars("data/processed_data/firEUrisk_wild_veg.tif") |> 
  get_large_polygons() 
terra::writeRaster(firEUrisk_wild_veg_large, 
                   "data/processed_data/firEUrisk_wild_veg_large.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Distances to large vegetated areas ----
set.values(firEUrisk_wild_veg_large)
firEUrisk_wild_veg_large_dist <- terra::distance(firEUrisk_wild_veg_large)
terra::writeRaster(firEUrisk_wild_veg_large_dist, 
                   "data/processed_data/firEUrisk_wild_veg_large_dist.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)

# Percentage-vegetation ----
firEUrisk_wv_perc <- get_perc_veg_buffer(
  .veg = firEUrisk_wild_veg,
  .radio_buffer = radio_buffer)
terra::writeRaster(firEUrisk_wv_perc, 
                   "data/processed_data/firEUrisk_wv_perc.tif",
                   filetype = "GTiff", 
                   overwrite = TRUE)
# Plot ----
levels(firEUrisk_fuel_250m) <- data.frame(id = firEUrisk_fuel_type$code,
                                          type =  firEUrisk_fuel_type$description)

plot(region_r, col = "lightgrey", axes = FALSE, legend = FALSE)
plot(firEUrisk_fuel_250m, axes = FALSE, add = T)
Figure 3.10: Wildland vegetation in Norway by type of vegetation (Source EFFIS Fuel Map)
Show the code 
# Read Effis data
firEUrisk_wild_veg <- terra::rast("data/processed_data/firEUrisk_wild_veg.tif")
dist_path <- "data/processed_data/firEUrisk_wild_veg_large_dist.tif"
firEUrisk_wild_veg_large_dist_2400 <- terra::rast(dist_path) |> 
  get_dist_max(2400)

# Plot
par(mfrow = c(1,2), adj = 0)

# Wildand vegetation (Yes = 1, No = NA)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(firEUrisk_wild_veg, col = "#009E73", axes = F, legend = F, add = T)
title("A)", adj = 0)
# Within a distance of 2.4 km (Yes = 1, No = 0)
plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(firEUrisk_wild_veg_large_dist_2400, col = c(NA, "#0072B2"), axes = F, legend = F, add = T)
title("B)", adj = 0)
Figure 3.11: a) Wildland vegetation in Norway, and b) grid cells within a distance to large vegetation areas (>5km2) of 2.4 km (Source EFFIS Fuel Map)

3.3 WUI classification

We have used the definition described by (Stewart et al. 2007; Bar-Massada 2021; Bar-Massada et al. 2013), which is shown in Figure 3.12:

Show the code 
knitr::include_graphics("figures/wui_map_approach_bar_massada_2021.png")
Figure 3.12: Zonal-based WUI mapping approach

We carried out a point-based approach with grid cells of 250 m x 250 m, and a buffer radio of 500 m. Our workflow was therefore:

flowchart TD
    A["Grid 250m x 250m"] --> B{"Building density > 6.17 #/km2"}
    B -- Yes --> C{"Vegetation > 50%"}
    B -- No --> D["No WUI"]
    C -- Yes --> E["Intermix"]
    C -- No --> F{"Within 2.4 km of large 
    forest (>5 km)"}
    F -- Yes --> G["Interface"]
    F -- No --> H["No WUI"]
    style D fill:lightgrey
    style E fill:#003366,color:#FFFFFF
    style G fill:#FF6633
    style H fill:lightgrey

Therefore, we did the following steps:

  1. Create a raster with the potential WUI cells (wui_potential)

  • Building density (in the buffer area) ≤ 6.16 #/km^2 –> NA

  • Building density (in the buffer area) > 6.16 #/km^2 –> 1

  1. Create a raster with the grids where the percentage of fuel vegetation in the buffer area is higher than 50% (wv_perc)
  • Percentage < 50% –> 0
  • Percentage ≥ 50% –> 1
  1. Create a raster indicating if the grid cell is within 2.4 km of a large (>5 km2) fuel vegetation area (wild_veg_large_dist_*)
  • Yes: 1
  • No: 0
  1. Sum both rasters (wv_perc + wui_potential), and generate a new raster (wui) with three options :

  • wui = NA –> No WUI

  • wui = 2 –> WUI intermix

  • wui = 1 –> Potential interface (it is reassigned to 0 in the raster)

  1. Sum this raster (wui) and wild_veg_large_dist_, we get therefore that:
  • wui = NA –> No WUI
  • wui = 0 –> No WUI
  • wui = 1 –> Interface
  • wui = 2 –> Intermix
  • wui = 3 –> Intermix
Show the code 
# Minimum number of buildings (N) in the buffer area ----
# WUI if building density (d) > 6.17 [#/km^2]
# N = d [#/km^2] * A [km^2]
# A = pi * radio^2
d_min <- 6.17
N_min <- round(d_min * area_buffer, 3)

wui_potential <- build_focal
wui_potential[wui_potential <= N_min] <- NA
wui_potential[wui_potential >  N_min] <- 1

plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(wui_potential, axes = F, legend = F, add = T)
Figure 3.13: Potential WUI grids based on building density
Show the code 
effis_wv_perc     <- terra::rast("data/processed_data/effis_wv_perc.tif")
corine_wv_perc    <- terra::rast("data/processed_data/corine_wv_perc.tif")
sr16_wv_perc      <- terra::rast("data/processed_data/sr16_wv_perc.tif")
nmd_wv_perc       <- terra::rast("data/processed_data/nmd_wv_perc.tif")
firEUrisk_wv_perc <- terra::rast("data/processed_data/firEUrisk_wv_perc.tif")
Show the code 
# EFFIS
effis_wui <- get_wui_raster(.veg_perc = effis_wv_perc,
                            .veg_dist = effis_wild_veg_large_dist_2400)

# CORINE
corine_wui <- get_wui_raster(.veg_perc = corine_wv_perc,
                             .veg_dist = corine_wild_veg_large_dist_2400)

# National maps
sr16_wui <- get_wui_raster(.veg_perc = sr16_wv_perc,
                           .veg_dist = sr16_wild_veg_large_dist_2400,
                           .region = no_v)

nrm_wui <- get_wui_raster(.veg_perc = nmd_wv_perc,
                          .veg_dist = nmd_wild_veg_large_dist_2400,
                          .region = se_v)

national_wui <- mosaic(sr16_wui, nrm_wui, fun = "max") |> 
  add_colours_wui()

# firEUrisk
firEUrisk_wui <- get_wui_raster(.veg_perc = firEUrisk_wv_perc,
                                .veg_dist = firEUrisk_wild_veg_large_dist_2400)

# Put all rasters in one 
wui_r <- c(effis_wui, corine_wui, national_wui, firEUrisk_wui)
names(wui_r) <- c("effis", "corine", "national", "firEUrisk")

terra::saveRDS(wui_r, "data/processed_data/wui_r.rds")

Map of WUI areas (Figure 3.14; Table 3.2).

Show the code 
wui_r <- terra::readRDS("data/processed_data/wui_r.rds")
Show the code 
f_plot <- function(.x){ plot(.x, 
                             axes = F, 
                             mar = c(0, 0, 0, 0),
                             legend = "topleft")  }

f_plot(wui_r$effis)
f_plot(wui_r$corine)
f_plot(wui_r$firEUrisk)
f_plot(wui_r$national)
(a) EFFIS
(b) CORINE
(c) FirEUrisk
(d) National databases
Figure 3.14: WUI maps of Norway based on different fuel/vegetation maps (resolution of 250 x 250 m)

Percentages at national scale from the different vegetation/fuel types (based on the number of grids cells of 250x 250 [m] for each WUI type):

Show the code 
# summary table 
tbl_summary_vegetation <- function(.x){
  
  .x |> 
    as_tibble() |> 
    drop_na() |> 
    pivot_longer(everything(), 
                 names_to = "name",
                 values_to = "value") |> 
    group_by(name, value) |> 
    summarise(count = n()) |> 
    ungroup() |> 
    group_by(name) |> 
    mutate(total = sum(count),
           pct = 100*count/total) |> 
    pivot_wider(names_from = "value", 
                values_from = c("count", "pct")) |> 
    mutate(pct_wui = pct_Interface + pct_Intermix) |> 
    select(starts_with("pct")) |> 
    relocate(pct_wui, .after = pct_No_wui) |> 
    gt(row_group_as_column = TRUE) |> 
    fmt_number(decimals = 1) |> 
    cols_label(
      name = "",
      pct_No_wui = md("**No WUI**"),
      pct_Interface = md("**Interface**"),
      pct_Intermix = md("**Intermix**"),
      pct_wui = md("**Total**")
    )  |> 
    tab_spanner(
      label = md("**WUI**"),
      columns = c(pct_Interface, pct_Intermix, pct_wui)
    ) |> 
    tab_options(
      table.font.size = 9,
      row_group.padding = px(1),
      row_group.font.weight = "bold",
      data_row.padding = px(1)
    )
  
}


# Summary WUI by country 
wui_r_no <- mask(wui_r, no_v)
wui_r_se <- mask(wui_r, se_v)

# Tables
tbl_summary_vegetation(wui_r_no) 
tbl_summary_vegetation(wui_r_se)
Table 3.2: Percentages of WUI by vegetation and country
(a) Norway

No WUI

WUI

Interface

Intermix

Total
corine 83.1 6.6 10.3 16.9
effis 82.5 5.9 11.6 17.5
firEUrisk 86.5 3.8 9.7 13.5
national 88.7 7.5 3.8 11.3
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
corine 91.2 4.1 4.7 8.8
effis 91.1 3.7 5.2 8.9
firEUrisk 92.7 2.2 5.1 7.3
national 90.6 4.5 4.9 9.4

Values by counties:

Show the code 
# Function for printing the table
tbl_summary_county <- function(.x) {
  .x |> 
    select(NUTS_NAME, name, starts_with("pct")) |> 
    gt(rowname_col = "name",
       groupname_col = "NUTS_NAME",
       row_group_as_column = FALSE) |> 
    fmt_number(decimals = 1) |> 
    cols_label(
      name = "",
      pct_No_wui = md("**No WUI**"),
      pct_Interface = md("**Interface**"),
      pct_Intermix = md("**Intermix**"),
      pct_wui = md("**Total**")
    )  |> 
    tab_spanner(
      label = md("**WUI**"),
      columns = c(pct_Interface, pct_Intermix, pct_wui)
    ) |> 
    tab_options(
      table.font.size = 9,
      row_group.background.color = "lightgrey",
      row_group.padding = px(1),
      row_group.font.weight = "bold",
      data_row.padding = px(1)
    ) 
}


# Summary by county
tbl_wui_county <- terra::extract(wui_r, counties, touches = TRUE) |>  
  pivot_longer(cols = -ID, 
               names_to = "name",
               values_to = "value") |> 
  group_by(ID, name, value) |> 
  summarise(count = n()) |> 
  ungroup() |> 
  group_by(ID, name) |> 
  mutate(total = sum(count), 
         pct = 100*count/total) |> 
  ungroup() |> 
  pivot_wider(names_from = "value",
              values_from = c("count", "pct")) |> 
  mutate(pct_wui = pct_Interface + pct_Intermix) 
  
# Add values to sf counties
counties_wui <- counties |> 
  st_as_sf()  %>% 
  mutate(ID := seq_len(nrow(.)))  |> 
  left_join(tbl_wui_county, by = "ID") |> 
  st_drop_geometry() |> 
  as_tibble() 

# Table by country & county 
filter(counties_wui, CNTR_CODE == "NO") |> 
  tbl_summary_county()
filter(counties_wui, CNTR_CODE == "SE") |> 
  tbl_summary_county()
Table 3.3: Percentages of WUI by county
(a) Norway

No WUI

WUI

Interface

Intermix

Total
Viken
corine 64.1 13.9 22.0 35.9
effis 64.6 12.4 22.9 35.4
firEUrisk 69.5 7.8 22.7 30.5
national 69.6 20.0 10.3 30.4
Rogaland
corine 81.2 9.9 8.9 18.8
effis 78.1 10.3 11.6 21.9
firEUrisk 87.2 5.5 7.4 12.8
national 91.2 6.4 2.4 8.8
Vestland
corine 83.5 6.9 9.6 16.5
effis 82.0 6.8 11.3 18.0
firEUrisk 88.3 4.1 7.7 11.7
national 93.0 4.6 2.4 7.0
Trøndelag
corine 83.6 8.4 8.0 16.4
effis 82.3 6.8 10.9 17.7
firEUrisk 87.2 4.7 8.2 12.8
national 90.9 7.0 2.2 9.1
Troms og Finnmark
corine 94.0 2.5 3.5 6.0
effis 93.9 2.4 3.8 6.1
firEUrisk 95.7 1.5 2.8 4.3
national 97.0 2.2 0.8 3.0
Innlandet
corine 78.3 7.1 14.6 21.7
effis 78.3 5.2 16.5 21.7
firEUrisk 80.5 4.1 15.4 19.5
national 82.3 11.4 6.3 17.7
Oslo
corine 66.0 12.5 21.5 34.0
effis 67.0 11.4 21.6 33.0
firEUrisk 70.9 7.3 21.8 29.1
national 64.6 23.8 11.6 35.4
Vestfold og Telemark
corine 74.9 7.5 17.6 25.1
effis 75.2 6.8 18.0 24.8
firEUrisk 78.6 4.4 17.0 21.4
national 79.6 12.1 8.3 20.4
Agder
corine 75.9 5.0 19.1 24.1
effis 76.0 4.5 19.5 24.0
firEUrisk 79.1 2.9 18.0 20.9
national 83.1 10.0 6.9 16.9
Møre og Romsdal
corine 78.9 9.6 11.5 21.1
effis 77.4 9.5 13.1 22.6
firEUrisk 85.3 5.9 8.8 14.7
national 84.8 10.8 4.5 15.2
Nordland
corine 88.4 5.2 6.4 11.6
effis 87.5 5.1 7.4 12.5
firEUrisk 92.5 2.7 4.9 7.5
national 94.8 3.7 1.5 5.2
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
Kalmar län
corine 88.5 4.3 7.2 11.5
effis 86.9 4.7 8.4 13.1
firEUrisk 89.6 2.4 8.0 10.4
national 87.7 4.8 7.6 12.3
Blekinge län
corine 86.8 5.6 7.6 13.2
effis 86.5 5.5 8.0 13.5
firEUrisk 90.3 2.7 7.1 9.7
national 86.3 6.5 7.2 13.7
Gävleborgs län
corine 90.8 3.7 5.6 9.2
effis 90.9 3.2 5.9 9.1
firEUrisk 92.0 2.4 5.7 8.0
national 90.4 3.9 5.8 9.6
Västerbottens län
corine 96.6 1.3 2.1 3.4
effis 96.7 1.1 2.3 3.3
firEUrisk 96.9 0.8 2.3 3.1
national 96.6 1.2 2.2 3.4
Norrbottens län
corine 98.3 0.8 1.0 1.7
effis 98.2 0.6 1.2 1.8
firEUrisk 98.5 0.4 1.1 1.5
national 98.2 0.8 1.0 1.8
Södermanlands län
corine 75.8 13.5 10.6 24.2
effis 75.8 12.6 11.6 24.2
firEUrisk 80.4 6.9 12.7 19.6
national 74.1 14.9 11.0 25.9
Värmlands län
corine 86.1 5.9 8.0 13.9
effis 86.1 5.4 8.5 13.9
firEUrisk 87.4 3.6 9.0 12.6
national 85.6 5.9 8.5 14.4
Uppsala län
corine 77.6 11.6 10.8 22.4
effis 77.4 10.9 11.7 22.6
firEUrisk 82.3 6.6 11.2 17.7
national 75.5 12.6 11.9 24.5
Kronobergs län
corine 88.2 3.6 8.2 11.8
effis 88.2 3.1 8.7 11.8
firEUrisk 89.0 2.1 8.9 11.0
national 87.9 3.8 8.2 12.1
Skåne län
corine 86.9 7.2 5.9 13.1
effis 84.3 8.1 7.6 15.7
firEUrisk 90.6 3.4 6.0 9.4
national 84.2 9.8 6.0 15.8
Stockholms län
corine 54.6 24.0 21.4 45.4
effis 54.2 21.9 23.8 45.8
firEUrisk 65.5 11.8 22.6 34.5
national 49.7 26.2 24.1 50.3
Västmanlands län
corine 87.2 8.0 4.8 12.8
effis 87.2 7.5 5.3 12.8
firEUrisk 90.5 4.5 5.1 9.5
national 85.6 9.2 5.2 14.4
Hallands län
corine 90.7 4.8 4.5 9.3
effis 90.4 4.8 4.9 9.6
firEUrisk 92.9 2.3 4.8 7.1
national 89.2 6.3 4.4 10.8
Västra Götalands län
corine 74.4 12.8 12.8 25.6
effis 74.3 11.9 13.8 25.7
firEUrisk 80.2 6.6 13.1 19.8
national 73.1 14.8 12.2 26.9
Jämtlands län
corine 97.7 0.9 1.4 2.3
effis 97.7 0.7 1.6 2.3
firEUrisk 97.9 0.6 1.6 2.1
national 97.6 1.0 1.4 2.4
Östergötlands län
corine 87.5 6.3 6.1 12.5
effis 87.1 6.2 6.8 12.9
firEUrisk 89.9 3.2 6.9 10.1
national 85.8 7.6 6.5 14.2
Västernorrlands län
corine 90.9 3.3 5.8 9.1
effis 91.0 3.0 6.0 9.0
firEUrisk 92.1 2.2 5.7 7.9
national 90.7 3.3 6.0 9.3
Dalarnas län
corine 92.9 3.0 4.1 7.1
effis 93.0 2.7 4.3 7.0
firEUrisk 93.9 1.7 4.4 6.1
national 92.4 3.1 4.4 7.6
Jönköpings län
corine 89.4 4.5 6.1 10.6
effis 89.5 3.8 6.7 10.5
firEUrisk 90.3 2.4 7.3 9.7
national 88.9 4.9 6.2 11.1
Örebro län
corine 81.8 7.7 10.5 18.2
effis 81.8 7.2 11.0 18.2
firEUrisk 85.5 3.9 10.6 14.5
national 80.2 9.0 10.8 19.8
Gotlands län
corine 82.1 11.0 6.9 17.9
effis 80.8 9.5 9.7 19.2
firEUrisk 85.4 6.0 8.6 14.6
national 81.3 12.2 6.4 18.7
Show the code 
# Add county ID
counties_cent_ID <- counties |>
  st_centroid() %>%
  mutate(x = st_coordinates(.)[,1],
         y = st_coordinates(.)[,2]) |> 
  arrange(desc(y), x) |> 
  mutate(ID = seq_along(counties$NUTS_NAME)) 

counties_ID <- counties |> 
  as_tibble() |> 
  left_join(counties_cent_ID, "NUTS_NAME") |> 
  st_sf()

# geofacet

mygrid <- data.frame(
  code = c("1 Troms og F.", 
           "2 Norrbotten", 
           "3 Nordland", 
           "4 Västerbotten", 
           "5 Trøndelag", 
           "7 Västernorrland", 
           "6 Jämtland", 
           "8 Møre og R.", 
           "10 Innlandet", 
           "11 Dalarna", 
           "9 Gävleborg", 
           "12 Vestland",
           "15 Oslo",
           "14 Viken",
           "13 Uppsala",
           "16 Västmanland",
           "20 Örebro", 
           "17 Värmland",
           "22 Rogaland", 
           "19 Vestfold og T.",
           "18 Stockholm", 
           "21 Södermanland",
           "24 Östergötland",
           "25 Västra Götaland",
           "23 Agder", 
           "29 Halland", 
           "27 Jönköping", 
           "28 Kalmar", 
           "26 Gotland", 
           "30 Kronoberg", 
           "31 Blekinge",
           "32 Skåne"),
  row = c(2, 2, 3, 3, 4, 4, 4, 5, 5, 5,
          5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 
          7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10),
  col = c(6, 7, 6, 7, 5, 7, 6, 4, 5, 6,
          7, 3, 5, 4, 9, 8, 7, 6, 3, 4, 
          9, 8, 7, 6, 4, 6, 7, 8, 9, 6, 7, 6),
  name = c("Troms og Finnmark",
           "Norrbotten", 
           "Nordland",
           "Västerbotten", 
           "Trøndelag", 
           "Västernorrland",
           "Jämtland", 
           "Møre og Romsdal",
           "Innlandet", 
           "Dalarna", 
           "Gävleborg",
           "Vestland", 
           "Oslo", 
           "Viken", 
           "Uppsala",
           "Västmanland",
           "Örebro", 
           "Värmland",
           "Rogaland", 
           "Vestfold og Telemark",
           "Stockholm", 
           "Södermanland", 
           "Östergötland", 
           "Västra Götaland", 
           "Agder", 
           "Halland", 
           "Jönköping", 
           "Kalmar",
           "Gotland", 
           "Kronoberg", 
           "Blekinge", 
           "Skåne"),
  stringsAsFactors = FALSE
)

# Maps

p_map <- ggplot() +
  geom_sf(data = counties_ID, 
          fill = "#FAFAFA",
          col = "lightgrey",
          size = 0.5) +
  geom_sf_text(data = counties_ID, 
               aes(label = ID),
               col = "red",
               size = 3) +
  theme_void()

p_bar <- counties_wui |> 
  select(NUTS_NAME, name, starts_with("pct")) |> 
  select(-pct_wui) |> 
  pivot_longer(cols = starts_with("pct"), 
               names_to = "type",
               values_to = "value") |> 
  mutate(type = str_remove_all(type, "pct_"),
         NUTS_NAME = str_remove_all(NUTS_NAME, "s län"),
         NUTS_NAME = str_remove_all(NUTS_NAME, " län")) |> 
  mutate(name = case_when(name == "corine" ~ "C",
                          name == "effis" ~ "E",
                          name == "national" ~ "N",
                          name == "firEUrisk" ~ "F")) |>
  mutate(name = factor(name, levels = c("N", "C", "E", "F"))) |> 
  mutate(type = factor(type, levels = c("No_wui", "Interface", "Intermix"))) |> 
  ggplot() +
  geom_bar(aes(x = name,
               y = value,
               fill = type),
           col = "grey",
           stat = "identity",
           position = "fill"
           ) +
  scale_fill_manual(values = c("#FAFAFA", "#FF6633", "#003366")) +
  geofacet::facet_geo(~NUTS_NAME,
                      grid = mygrid, 
                      label = "code"
                      ) +
  labs(x = "",
       y = "") +
  theme_bw() +
  theme(
    strip.text = element_text(size = 9, face = "bold")
  )

ggsave(plot = p_map, 
       filename = "figures/counties_map.png", 
       width = 3, height = 3)

ggsave(plot = p_bar, 
       filename = "figures/counties_wui_perc.png", 
       width = 15, height = 20)

WUI percentage by county

Area of wild vegetation by county (Table 3.4): we calculate the number of vegetation cells (i.e., grid centroid is inside the polygon), and calculate the area based on the raster resolution.

Show the code 
# Number of buildings by county  ----
county_build_count <- build_counts |> 
  # Number of buildings per county
  extract(y = vect(counties),
          fun = "sum",
          na.rm = TRUE,
          touches = TRUE,
          bind = TRUE) |> 
  # sf object
  st_as_sf() %>% 
  # Area polygons
  mutate(area = units::set_units(st_area(.), "km2"))  |> 
  # Building density by county
  mutate(build_den = sum/area)

# Percentage vegetation by county ----

# Function for getting vegetation area by counties 
f_pct_veg <- function(df) {
  df |> # sf object with vegetation data  
    # Number of cells
    extract(y = vect(counties),
            fun = "sum",
            na.rm = TRUE,
            touches = TRUE,
            bind = TRUE) |>  
    st_as_sf() |> 
    # Area (number cells x area cells)
    mutate(
      area_km = wild_veg * ((res(df))[1]/1000 * (res(df ))[2]/1000)
    ) |> 
    # Get only values 
    st_drop_geometry() |> 
    select(NUTS_NAME, CNTR_CODE, area_km)
}


## Run for each vegetation type

# Effis
effis_county  <- f_pct_veg(effis_wild_veg) |> 
    rename(area_km_effis = area_km)

# Corine
corine_county <- f_pct_veg(corine_wild_veg) |> 
  rename(area_km_corine = area_km)   

# National
sr16_county <- f_pct_veg(sr16_wild_veg) |> 
  filter(CNTR_CODE  == "NO")
nmd_county  <- f_pct_veg(nmd_wild_veg)  |> 
  filter(CNTR_CODE  == "SE")

national_county <- rbind(sr16_county, nmd_county) |> 
  rename(area_km_national = area_km)

# FirEUrisk
firEUrisk_county <- f_pct_veg(firEUrisk_wild_veg) |> 
  rename(area_km_firEUrisk = area_km)   

# Joint tables
county_veg <- national_county |> 
  left_join(corine_county) |> 
  left_join(effis_county) |> 
  left_join(firEUrisk_county)

# Add vegetation to the summary table ----
df_summary <- county_build_count |> 
  st_drop_geometry() |> 
  left_join(county_veg) |> 
  mutate(across(.cols = starts_with("area_km"),
                .fns = ~ 100 * (.x / as.numeric(area)),
                .names = "perc_{.col}"))

# Summary tables by country and county ----
summary_country <- df_summary |> 
  select(CNTR_CODE, sum, area, starts_with("area_km")) |> 
  group_by(CNTR_CODE) |> 
  summarise_all(sum) |> 
  mutate(NUTS_NAME = "Total",
         build_den = sum/area) |> 
  mutate(across(.cols = starts_with("area_km"),
                .fns = ~ 100 * (.x / as.numeric(area)),
                .names = "perc_{.col}")) |> 
  relocate(build_den, .after = area) |> 
  relocate(NUTS_NAME, .after = CNTR_CODE) |> 
  select(CNTR_CODE, NUTS_NAME, area,
         sum, build_den,
         starts_with("perc")) |> 
  mutate(area = as.numeric(area),
         build_den = as.numeric(build_den),
         CNTR_CODE = case_when(CNTR_CODE == "NO" ~ "NORWAY",
                               CNTR_CODE == "SE" ~ "SWEDEN"))
  
summary_county <- df_summary |> 
  select(CNTR_CODE, NUTS_NAME, area,
         sum, build_den,
         starts_with("perc")) |>
  arrange(CNTR_CODE, NUTS_NAME) |> 
  mutate(area = as.numeric(area),
         build_den = as.numeric(build_den),
         CNTR_CODE = case_when(CNTR_CODE == "NO" ~ "NORWAY",
                               CNTR_CODE == "SE" ~ "SWEDEN"))

# Print table ----
summary_county |> 
  rbind(summary_country) |> 
  gt(groupname_col = "CNTR_CODE",
     row_group_as_column = FALSE) |> 
  fmt_number(c(build_den, starts_with("perc")), decimals = 1) |> 
  fmt_number(c(sum, area), decimals = 0) |> 
  tab_options(
      table.font.size = 9,
      row_group.background.color = "lightgrey",
      row_group.padding = px(1),
      row_group.font.weight = "bold",
      data_row.padding = px(1)
    ) |> 
  tab_spanner(
      label = md("**Area covered by vegetation [%]**"),
      columns = starts_with("perc")
    ) |> 
  tab_spanner(
      label = md("**Buildings**"),
      columns = c(sum, build_den)
    ) |> 
  cols_label(
      NUTS_NAME = md("**County**"),
      area = md("**Area [km2]**"),
      sum = md("**Count**"),
      build_den = md("**Density [#/km2]**"),
      perc_area_km_national = md("**National**"),
      perc_area_km_effis = md("**Effis**"),
      perc_area_km_corine = md("**Corine**"),
      perc_area_km_firEUrisk = md("**FirEurisk**")
    ) |> 
  tab_style(
    style = cell_text(weight = "bold"),  # Apply bold text
    locations = cells_body(
      rows = c(33, 34)  # specify the row number
    )
  )
Table 3.4: Number of buildings and wildland vegetation percentage by county (and total)

County

Area [km2]

Buildings

Area covered by vegetation [%]

Count

Density [#/km2]

National

Corine

Effis

FirEurisk
NORWAY
Agder 16,428 272,049 16.6 31.5 73.5 89.2 70.7
Innlandet 52,068 579,653 11.1 31.7 54.7 83.8 71.3
Møre og Romsdal 14,973 234,325 15.6 19.3 48.7 74.5 55.2
Nordland 37,972 254,941 6.7 14.1 40.2 75.7 62.8
Oslo 455 129,178 284.0 42.5 62.7 64.3 67.7
Rogaland 9,352 303,296 32.4 16.1 43.3 73.7 49.6
Troms og Finnmark 74,774 232,001 3.1 9.6 40.7 83.6 67.5
Trøndelag 41,480 410,427 9.9 17.8 49.6 86.3 52.0
Vestfold og Telemark 17,449 364,558 20.9 37.2 72.1 85.2 68.3
Vestland 33,808 475,898 14.1 15.9 42.0 73.5 58.6
Viken 24,575 869,036 35.4 38.2 66.8 77.5 68.7
Total 323,335 4,125,362 12.8 20.9 50.0 81.2 63.9
SWEDEN
Blekinge län 3,035 41,010 13.5 71.6 77.8 80.8 79.2
Dalarnas län 30,395 133,414 4.4 83.5 81.5 89.4 81.7
Gotlands län 3,179 29,900 9.4 49.5 51.5 60.5 52.5
Gävleborgs län 19,718 100,559 5.1 84.9 83.8 86.9 86.5
Hallands län 5,713 72,407 12.7 61.8 61.9 65.7 66.8
Jämtlands län 54,089 61,692 1.1 72.2 79.7 90.0 69.1
Jönköpings län 11,745 77,384 6.6 71.6 72.0 76.6 81.9
Kalmar län 11,632 96,796 8.3 72.9 74.4 79.9 80.5
Kronobergs län 9,430 59,661 6.3 76.1 76.3 81.0 84.9
Norrbottens län 105,907 105,850 1.0 62.7 71.0 89.4 64.7
Skåne län 11,363 420,689 37.0 37.8 36.9 42.9 37.1
Stockholms län 7,073 440,098 62.2 59.0 55.8 59.1 54.0
Södermanlands län 7,047 76,316 10.8 56.1 55.2 57.5 60.7
Uppsala län 8,618 107,607 12.5 66.0 63.6 66.8 64.0
Värmlands län 21,910 144,947 6.6 71.4 69.7 73.0 74.4
Västerbottens län 59,240 107,866 1.8 78.4 80.8 91.0 76.5
Västernorrlands län 23,070 127,289 5.5 88.0 85.9 89.9 88.7
Västmanlands län 5,689 50,536 8.9 64.0 59.3 64.5 62.4
Västra Götalands län 28,859 542,657 18.8 53.4 54.0 57.4 57.1
Örebro län 9,687 103,333 10.7 70.2 69.5 72.6 72.7
Östergötlands län 12,257 137,746 11.2 59.0 59.2 61.6 64.5
Total 449,657 3,037,757 6.8 69.5 72.3 81.4 70.6

Total area of WUI (by type of vegetation):

Show the code 
wui_r |>  
  as_tibble() |> 
  drop_na() |> 
  pivot_longer(col = everything(),
               names_to = "name") |> 
  group_by(name, value) |>
  summarise(n = n()) |> 
  ungroup() |> 
  mutate(area = n * (0.25 * 0.25)) |> #km2
  pivot_wider(names_from = name, 
              values_from = c(n, area)) |> 
  select(value, starts_with("area")) |> 
  gt() |>
  fmt_number(decimals = 0) |> 
  cols_label(
    value = "",
    area_corine = md("**corine**"),
    area_effis = md("**effis**"),
    area_firEUrisk ~ md("**firEUrisk**"),
    area_national ~ md("**national**")
  ) |> 
  tab_options(
    table.font.size = 9,
    row_group.background.color = "lightgrey",
    row_group.padding = px(1),
    row_group.font.weight = "bold",
    data_row.padding = px(1)
  )

corine

effis

firEUrisk

national
No_wui 686,335 683,812 704,526 702,214
Interface 40,111 36,302 22,444 45,065
Intermix 55,420 61,752 54,896 34,588

Number of buildings in WUI

Show the code 
c(build_counts, wui_r) |> 
  as_tibble() |>   
  mutate(sum = ifelse(is.na(sum), 0, sum)) |> 
  drop_na() |> 
  pivot_longer(col = -sum, names_to = "name") |> 
  group_by(name, value) |>
  summarise(buildings = sum(sum)) |> 
  ungroup() |> 
  pivot_wider(names_from = name, 
              values_from = buildings) |> 
  gt(rowname_col =  "value") |> 
  grand_summary_rows(
    fns = list(Total = ~sum(.))
  ) |> 
   cols_label(
    corine = md("**corine**"),
    effis = md("**effis**"),
    firEUrisk ~ md("**firEUrisk**"),
    national ~ md("**national**")
  ) |>
  tab_options(
    table.font.size = 9,
    row_group.padding = px(1),
    row_group.font.weight = "bold",
    data_row.padding = px(1)
  )
Table 3.5: Number of buildings in each WUI classification

corine

effis

firEUrisk

national
No_wui 2272345 2185435 3662215 2915152
Interface 3449156 3213257 1579839 3460202
Intermix 1432454 1755263 1911901 778601
Total 7153955 7153955 7153955 7153955

3.3.1 Effect distance to large forest

Show the code 
# Distance to large forest areas 

# Distance variation
max_distances <- c(200, 500, 1000, 2000, 2400)

# WUIs
f_wui <- function(.x, .y){
  
  veg_perc <- .x
  
  get_wui_raster(.veg_perc = veg_perc,
                 .veg_dist = .y,
                 .pot_wui = wui_potential,
                 .region = region_v)
  
}

# Summary table
tbl_summary_distance <- function(.x){

.x |>
  as_tibble() |>
  drop_na() |>
  pivot_longer(everything(),
               names_to = "name",
               values_to = "value") |>
  group_by(name, value) |>
  summarise(count = n()) |>
  ungroup() |>
  group_by(name) |>
  mutate(total = sum(count),
         pct = 100*count/total) |>
  pivot_wider(names_from = "value",
              values_from = c("count", "pct")) |>
  mutate(pct_wui = pct_Interface + pct_Intermix) |>
  select(starts_with("pct")) |>
  relocate(pct_wui, .after = pct_No_wui) |>
  mutate(name = factor(name,
                       levels = c("d200",
                                  "d500",
                                  "d1000",
                                  "d2000",
                                  "d2400"))) |>
  arrange(name) |>
  gt(row_group_as_column = TRUE) |>
  fmt_number(decimals = 1) |>
  cols_label(
    name = "",
    pct_No_wui = md("**No WUI**"),
    pct_Interface = md("**Interface**"),
    pct_Intermix = md("**Intermix**"),
    pct_wui = md("**Total**")
  )  |>
  tab_spanner(
    label = md("**WUI**"),
    columns = c(pct_Interface, pct_Intermix, pct_wui)
  ) |>
  tab_options(
    table.font.size = 9,
    row_group.padding = px(1),
    row_group.font.weight = "bold",
    data_row.padding = px(1)
  )

}

3.3.1.1 Effis

Show the code 
# Raster with distances to large forest (in m)
dist_path <- "data/processed_data/effis_wild_veg_large_dist.tif"
effis_wild_veg_large_dist <- rast(dist_path)

# Withing the max. distance
effis_dist_variation<- map2(effis_wild_veg_large_dist, max_distances, get_dist_max)
names(effis_dist_variation) <- c("d200", "d500", "d1000", "d2000", "d2400")

# Wui vs distance
wui_effis <- map2(effis_wv_perc,
                  effis_dist_variation,
                  f_wui)
wui_effis <- rast(wui_effis)
names(wui_effis) <- c("d200", "d500", "d1000", "d2000", "d2400")
set.values(wui_effis)
terra::saveRDS(wui_effis, "data/processed_data/wui_effis.rsd")

# Plot
panel(wui_effis, axes = FALSE)
Figure 3.15: WUI change based on EFFIS fuel map
Show the code 
wui_effis <- terra::readRDS("data/processed_data/wui_effis.rsd")

# Data by country
wui_effis_no <- mask(wui_effis, no_v)
wui_effis_se <- mask(wui_effis, se_v)
Show the code 
tbl_summary_distance(wui_effis_no)
tbl_summary_distance(wui_effis_se)
Table 3.6: Changes in WUI percentages of WUI by country (EFFIS)
(a) Norway

No WUI

WUI

Interface

Intermix

Total
d200 85.9 2.5 11.6 14.1
d500 83.6 4.7 11.6 16.4
d1000 83.0 5.4 11.6 17.0
d2000 82.6 5.8 11.6 17.4
d2400 82.5 5.9 11.6 17.5
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
d200 93.5 1.4 5.2 6.5
d500 92.1 2.7 5.2 7.9
d1000 91.6 3.3 5.2 8.4
d2000 91.2 3.6 5.2 8.8
d2400 91.1 3.7 5.2 8.9

3.3.1.2 Corine

Show the code 
# Raster with distances to large forest (in m)
dist_path <- "data/processed_data/corine_wild_veg_large_dist.tif"
corine_wild_veg_large_dist <- rast(dist_path)

# Withing the max. distance
corine_dist_variation<- map2(corine_wild_veg_large_dist, max_distances, get_dist_max)
names(corine_dist_variation) <- c("d200", "d500", "d1000", "d2000", "d2400")

# Wui vs distance
wui_corine <- map2(corine_wv_perc,
                   corine_dist_variation,
                   f_wui)
wui_corine <- rast(wui_corine)
names(wui_corine) <- c("d200", "d500", "d1000", "d2000", "d2400")
terra::saveRDS(wui_corine, "data/processed_data/wui_corine.rsd")

# Plot
panel(wui_corine, axes = FALSE)
Figure 3.16: WUI change based on CORINE vegetation map
Show the code 
wui_corine <- terra::readRDS("data/processed_data/wui_corine.rsd")

# Data by country
wui_corine_no <- mask(wui_corine, no_v)
wui_corine_se <- mask(wui_corine, se_v)
Show the code 
tbl_summary_distance(wui_corine_no)
tbl_summary_distance(wui_corine_se)
Table 3.7: Changes in WUI percentages of WUI by country (CORINE)
(a) Norway

No WUI

WUI

Interface

Intermix

Total
d200 86.9 2.8 10.3 13.1
d500 84.7 5.0 10.3 15.3
d1000 83.8 5.9 10.3 16.2
d2000 83.2 6.5 10.3 16.8
d2400 83.1 6.6 10.3 16.9
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
d200 93.7 1.6 4.7 6.3
d500 92.4 2.9 4.7 7.6
d1000 91.7 3.6 4.7 8.3
d2000 91.3 4.0 4.7 8.7
d2400 91.2 4.1 4.7 8.8

3.3.1.3 National datasets

Show the code 
# Norway ----
dist_path <- "data/processed_data/sr16_wild_veg_large_dist.tif"
sr16_wild_veg_large_dist <- rast(dist_path)
sr16_dist_variation<- map2(sr16_wild_veg_large_dist, max_distances, get_dist_max)
names(sr16_dist_variation) <- c("d200", "d500", "d1000", "d2000", "d2400")
wui_sr16 <- map2(sr16_wv_perc, sr16_dist_variation, f_wui)
wui_sr16 <- rast(wui_sr16)
names(wui_sr16) <- c("d200", "d500", "d1000", "d2000", "d2400")
terra::saveRDS(wui_sr16, "data/processed_data/wui_sr16.rsd")

# Sweden ----
dist_path <- "data/processed_data/nmd_wild_veg_large_dist.tif"
nmd_wild_veg_large_dist <- rast(dist_path)
nmd_dist_variation<- map2(nmd_wild_veg_large_dist, max_distances, get_dist_max)
names(nmd_dist_variation) <- c("d200", "d500", "d1000", "d2000", "d2400")
wui_nmd <- map2(nmd_wv_perc, nmd_dist_variation, f_wui)
wui_nmd <- rast(wui_nmd)
names(wui_nmd) <- c("d200", "d500", "d1000", "d2000", "d2400")
terra::saveRDS(wui_nmd, "data/processed_data/wui_nmd.rsd")

# Merge datasets ----
wui_national <- mosaic(wui_sr16, wui_nmd, fun = "max") |> 
  add_colours_wui()

# Add names
names(wui_national) <- c("d200", "d500", "d1000", "d2000", "d2400")
set.values(wui_national)
# Export
terra::saveRDS(wui_national, "data/processed_data/wui_national.rsd")

# Plot region  
panel(wui_national, axes = FALSE)
Figure 3.17: WUI change based on national vegetation maps
Show the code 
wui_national <- terra::readRDS("data/processed_data/wui_national.rsd")

# Data by country
wui_national_no <- mask(wui_national, no_v)
wui_national_se <- mask(wui_national, se_v)
Show the code 
tbl_summary_distance(wui_national_no)
tbl_summary_distance(wui_national_se)
Table 3.8: Changes in WUI percentages of WUI by country (national datasets - SR16 & NMD)
(a) Norway

No WUI

WUI

Interface

Intermix

Total
d200 94.6 1.6 3.8 5.4
d500 92.7 3.5 3.8 7.3
d1000 91.0 5.2 3.8 9.0
d2000 89.2 7.0 3.8 10.8
d2400 88.7 7.5 3.8 11.3
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
d200 93.7 1.5 4.9 6.3
d500 92.4 2.7 4.9 7.6
d1000 91.6 3.6 4.9 8.4
d2000 90.8 4.3 4.9 9.2
d2400 90.6 4.5 4.9 9.4

3.3.1.4 FirEUrisk

Show the code 
# Raster with distances to large forest (in m)
dist_path <- "data/processed_data/firEUrisk_wild_veg_large_dist.tif"
firEUrisk_wild_veg_large_dist <- rast(dist_path)

# Withing the max. distance
firEUrisk_dist_variation<- map2(firEUrisk_wild_veg_large_dist, max_distances, get_dist_max)
names(firEUrisk_dist_variation) <- c("d200", "d500", "d1000", "d2000", "d2400")

# Wui vs distance
wui_firEUrisk <- map2(firEUrisk_wv_perc,
                   firEUrisk_dist_variation,
                   f_wui)
wui_firEUrisk <- rast(wui_firEUrisk)
names(wui_firEUrisk) <- c("d200", "d500", "d1000", "d2000", "d2400")
terra::saveRDS(wui_firEUrisk, "data/processed_data/wui_firEUrisk.rsd")

# Plot
panel(wui_firEUrisk, axes = FALSE)
Figure 3.18: WUI change based on firEUrisk fuel map
Show the code 
wui_firEUrisk <- terra::readRDS("data/processed_data/wui_firEUrisk.rsd")

# Data by country
wui_firEUrisk_no <- mask(wui_firEUrisk, no_v)
wui_firEUrisk_se <- mask(wui_firEUrisk, se_v)
Show the code 
tbl_summary_distance(wui_firEUrisk_no)
tbl_summary_distance(wui_firEUrisk_se)
Table 3.9: Changes in WUI percentages of WUI by country (firEUrisk)
(a) Norway

No WUI

WUI

Interface

Intermix

Total
d200 88.4 1.9 9.7 11.6
d500 86.9 3.4 9.7 13.1
d1000 86.8 3.5 9.7 13.2
d2000 86.6 3.7 9.7 13.4
d2400 86.5 3.8 9.7 13.5
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
d200 93.8 1.1 5.1 6.2
d500 92.9 2.1 5.1 7.1
d1000 92.8 2.1 5.1 7.2
d2000 92.8 2.2 5.1 7.2
d2400 92.7 2.2 5.1 7.3

Aggregate all maps in one figure (Note: data are re-sampled to 500859 cells in the plots)

Show the code 
# Plot
p_effis <- ggplot() + 
  geom_spatraster(data = wui_effis, maxcell = 2000000) +
  facet_wrap(~lyr, ncol = 1) +
  labs(title = "Effis") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))
  
p_corine <- ggplot() + 
  geom_spatraster(data = wui_corine, maxcell = 2000000) +
  facet_wrap(~lyr, ncol = 1) +
  labs(title = "Corine") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

p_national <- ggplot() + 
  geom_spatraster(data = wui_national, maxcell = 2000000) +
  facet_wrap(~lyr, ncol = 1) +
  labs(title = "National") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

p_firEUrisk <- ggplot() + 
  geom_spatraster(data = wui_firEUrisk, maxcell = 2000000) +
  facet_wrap(~lyr, ncol = 1) +
  labs(title = "firEUrisk") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

p_effis + p_corine + p_national + p_firEUrisk +
  plot_layout(guides = "collect",
              ncol = 4)
Figure 3.19: WUI by vegetation type map and threshold distance

Plot variation in the percentage of WUI with the distance to large forest areas by each vegetation/fuel map.

Show the code 
get_tbl_pct <- function(.x) {
  
  t <- .x |>
    as_tibble() |> 
    drop_na() |> 
    pivot_longer(everything(),
                 names_to = "distance",
                 values_to = "value") |>
    group_by(distance, value) |>
    summarise(count = n()) |>
    ungroup() |> 
    group_by(distance) |>
    mutate(total = sum(count),
           pct = 100*count/total) |> 
    ungroup() |> 
    pivot_wider(names_from = "value",
                values_from = c("count", "pct")) |> 
    select(distance, starts_with("pct")) |> 
    mutate(distance = gsub("d", "", distance),
           distance = as.numeric(distance)) |> 
    arrange(distance)
  
  return(t)
  
}

# Combine rasters
raster_list <- list(wui_effis_no,
                    wui_corine_no,
                    wui_national_no,
                    wui_firEUrisk_no,
                    wui_effis_se,
                    wui_corine_se,
                    wui_national_se,
                    wui_firEUrisk_se)   

# Tables with percentages
dist_list <- raster_list |> purrr::map(get_tbl_pct)

# Add names to de list
names(dist_list) <- c("Norway_effis",
                      "Norway_corine",
                      "Norway_national",
                      "Norway_firEUrisk",
                      "Sweden_effis",
                      "Sweden_corine",
                      "Sweden_national",
                      "Sweden_firEUrisk")

# Merge with names
dist_tbl <- bind_rows(dist_list, .id = "Ticker") |> 
  separate(Ticker, c("country", "vegetation"))

# add value for 2400 m 
get_tbl_pct_2 <- function(.x){
  
t <- .x |>
  as_tibble() |> 
  drop_na() |> 
  pivot_longer(everything(),
               names_to = "vegetation",
               values_to = "value") |>
  group_by(vegetation, value) |>
  summarise(count = n()) |>
  ungroup() |> 
  group_by(vegetation) |>
  mutate(total = sum(count),
         pct = 100*count/total) |> 
  ungroup() |> 
  pivot_wider(names_from = "value",
              values_from = c("count", "pct")) |> 
  select(vegetation, starts_with("pct"))

return(t)

}
  
# Plot
ggplot(data = dist_tbl,
       aes(x = distance,
           y = pct_Interface,
           col = vegetation)) +
  geom_point() +
  geom_line() +
  facet_grid(~country) +
  labs(x = "Distance [m]",
       y = "Perc. [%]"
       # title = "Percentage of WUI Interface vs. distance to forest"
       ) +
  theme_bw() +
  theme(strip.text.x = element_text(size = 12, face = "italic"))
Figure 3.20: Change in the percentage of WUI with the theshold distance to large forest areas

3.3.2 Effect threshold building density

In the following figure we can see the distribution of buildings density (#/\(km^2\)). The number of buildings in the buffer area is then calculated multiplying the building density by the area of the buffer (i.e., minimum number of buildings (Nmin) = density threshold [#/\(km^2\)] * Buffer area [\(km^2\)]). In red is the standard threshold for classifying an area as possible WUI, then we have used half of this value (i.e., 3.08 #/\(km^2\)), representing approximately 25% of the data. Finally, we have also tested with 1 #/\(km^2\), which is equivalent to classified any area with buildings as potential WUI.

Show the code 
build_focal_no <- mask(build_focal, no_v) 

build_focal_se <- mask(build_focal, se_v)

# Get density table

get_d_counts <- function(.x) {
  
  build_focal |> 
    mask(.x) |> 
    as_tibble() |> 
    drop_na() |> 
    mutate(d_km2 = focal_sum / area_buffer)
  
}

df_counts <- list(no_v, se_v) |> 
  map(get_d_counts) 
names(df_counts) <- c("Norway", "Sweden")

df_counts_bind <- bind_rows(df_counts, .id = "country") 
  

# Create tables with summary statistics to add to the plot
get_tbl_quantiles <- function(.x) {
  
  tbl <- .x |> 
    select(d_km2) |>
    quantile(prob = seq(0, 1, 0.25), na.rm = T) |>
    round(2) |> 
    as_tibble() |> 
    rename(d = value) |> 
    mutate(Quantile = c("Min", "Q1", "Median", "Q3", "Max"),
           Quantile = forcats::fct_reorder(Quantile, d)) |> 
    select(Quantile, d)
  
  tbl_quantiles <- tibble(x = 2800,
                          y = 4.0e+5, 
                          table = list(tbl))
  
  return(tbl_quantiles)
  
}

tbl_quantiles <- df_counts |>
  map(get_tbl_quantiles) |> 
  bind_rows(.id = "country") 


# Plot
ggplot(df_counts_bind, aes(x = d_km2)) + 
  geom_histogram(bins = 50, color = "darkblue", fill = "lightblue") +
  facet_grid(~country) +
  scale_x_log10() +
  geom_vline(xintercept = 6.172,
             linetype = "dashed",
             color = "red",
             linewidth = 1.3) +
  annotate("text", x = 8, y = 420000,
           label = "6.172", colour = "red", parse = TRUE) + 
  geom_vline(xintercept = 3.08,
             linetype = "dashed",
             color = "orange",
             linewidth = 1.3) +
  annotate("text", x = 4, y = 420000,
           label = "3.08", colour = "orange", parse = TRUE) + 
  geom_vline(xintercept = 1.00,
             linetype = "dashed", 
             color = "blue",
             linewidth = 1.3) +
  annotate("text", x = 1.1, y = 420000,
           label = "1", colour = "blue", parse = TRUE) + 
  labs(title = "Building density",
       x = "#/km2",
       y = "") +
  theme_bw() +
  theme(legend.key = element_rect(fill = NA),
        legend.background = element_rect(color = "gray", linetype = "solid"),
        legend.position = c(.95, .95),
        legend.justification = c("right", "top"),
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6),
        strip.text.x = element_text(size = 12, face = "italic")) +
  # Add tabble with summary statistics
  ggpp::geom_table(data = tbl_quantiles,
                   aes(x = x, y = y, label = table),
                   table.theme = ggpp::ttheme_gtplain(base_size = 12))

Potential WUI based on different building density:

3.3.2.1 Building density 3.08 #/km2

Show the code 
# Minimum number of buildings (N) in the buffer area ----
# WUI if building density (d) > 6.17 [#/km^2]
# N = d [#/km^2] * A [km^2]
# A = pi * radio^2

N_min_3 <- round(3.08 * area_buffer, 3)

wui_potential_3 <- build_focal
wui_potential_3[wui_potential_3 <= N_min_3] <- NA
wui_potential_3[wui_potential_3 >  N_min_3] <- 1

plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(wui_potential_3, axes = F, legend = F, add = T)

Show the code 
# EFFIS
effis_wui_3 <- get_wui_raster(.veg_perc = effis_wv_perc,
                              .veg_dist = effis_wild_veg_large_dist_2400,
                              .pot_wui = wui_potential_3)

# CORINE
corine_wui_3 <- get_wui_raster(.veg_perc = corine_wv_perc,
                               .veg_dist = corine_wild_veg_large_dist_2400,
                               .pot_wui = wui_potential_3)

# National maps
sr16_wui_3 <- get_wui_raster(.veg_perc = sr16_wv_perc,
                             .veg_dist = sr16_wild_veg_large_dist_2400,
                             .pot_wui = wui_potential_3,
                             .region = no_v)

nrm_wui_3 <- get_wui_raster(.veg_perc = nmd_wv_perc,
                            .veg_dist = nmd_wild_veg_large_dist_2400,
                            .pot_wui = wui_potential_3,
                            .region = se_v)

national_wui_3 <- mosaic(sr16_wui_3, nrm_wui_3, fun = "max") |> 
  add_colours_wui()

# firEUrisk
firEUrisk_wui_3 <- get_wui_raster(.veg_perc = firEUrisk_wv_perc,
                                  .veg_dist = firEUrisk_wild_veg_large_dist_2400,
                                  .pot_wui = wui_potential_3)

# Put all rasters in one 
wui_r_3 <- c(effis_wui_3, corine_wui_3, national_wui_3, firEUrisk_wui_3)
names(wui_r_3) <- c("effis", "corine", "national", "firEUrisk")

terra::saveRDS(wui_r_3, "data/processed_data/wui_r_3.rds")

Map of WUI areas (Figure 3.21; Table 3.10).

Show the code 
wui_r_3 <- terra::readRDS("data/processed_data/wui_r_3.rds")
Show the code 
f_plot <- function(.x){ plot(.x, 
                             axes = F, 
                             mar = c(0, 0, 0, 0),
                             legend = "topleft")  }

f_plot(wui_r_3$effis)
f_plot(wui_r_3$corine)
f_plot(wui_r_3$firEUrisk)
f_plot(wui_r_3$national)
(a) EFFIS
(b) CORINE
(c) FirEUrisk
(d) National databases
Figure 3.21: WUI maps of Norway based on different fuel/vegetation maps (resolution of 250 x 250 m)
Show the code 
# Summary WUI by country
wui_r_3_no <- mask(wui_r_3, no_v)
wui_r_3_se <- mask(wui_r_3, se_v)

# summary table 
tbl_summary_vegetation(wui_r_3_no)
tbl_summary_vegetation(wui_r_3_se)
Table 3.10: Percentages of WUI by vegetation
(a) Norway

No WUI

WUI

Interface

Intermix

Total
corine 79.7 7.3 13.0 20.3
effis 79.0 6.3 14.8 21.0
firEUrisk 83.7 4.4 12.0 16.3
national 86.5 8.4 5.1 13.5
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
corine 89.3 4.5 6.1 10.7
effis 89.2 4.2 6.7 10.8
firEUrisk 91.0 2.5 6.5 9.0
national 88.7 5.0 6.3 11.3

3.3.2.2 Building density 1 #/km2

Equivalent to select any area with building as potential WUI.

Show the code 
# Minimum number of buildings (N) in the buffer area ----
# WUI if building density (d) > 6.17 [#/km^2]
# N = d [#/km^2] * A [km^2]
# A = pi * radio^2

N_min_1 <- round(1.00 * area_buffer, 3)

wui_potential_1 <- build_focal
wui_potential_1[wui_potential_1 <= N_min_1] <- NA
wui_potential_1[wui_potential_1 >  N_min_1] <- 1

plot(region_r, col = "lightgrey", axes = F, legend = F)
plot(wui_potential_1, axes = F, legend = F, add = T)

Show the code 
# EFFIS
effis_wui_1 <- get_wui_raster(.veg_perc = effis_wv_perc,
                              .veg_dist = effis_wild_veg_large_dist_2400,
                              .pot_wui = wui_potential_1)

# CORINE
corine_wui_1 <- get_wui_raster(.veg_perc = corine_wv_perc,
                               .veg_dist = corine_wild_veg_large_dist_2400,
                               .pot_wui = wui_potential_1)

# National maps
sr16_wui_1 <- get_wui_raster(.veg_perc = sr16_wv_perc,
                             .veg_dist = sr16_wild_veg_large_dist_2400,
                             .pot_wui = wui_potential_1,
                             .region = no_v)

nrm_wui_1 <- get_wui_raster(.veg_perc = nmd_wv_perc,
                            .veg_dist = nmd_wild_veg_large_dist_2400,
                            .pot_wui = wui_potential_1,
                            .region = se_v)

national_wui_1 <- mosaic(sr16_wui_1, nrm_wui_1, fun = "max") |> 
  add_colours_wui()

# firEUrisk
firEUrisk_wui_1 <- get_wui_raster(.veg_perc = firEUrisk_wv_perc,
                                  .veg_dist = firEUrisk_wild_veg_large_dist_2400,
                                  .pot_wui = wui_potential_1)

# Put all rasters in one 
wui_r_1 <- c(effis_wui_1, corine_wui_1, national_wui_1, firEUrisk_wui_1)
names(wui_r_1) <- c("effis", "corine", "national", "firEUrisk")

terra::saveRDS(wui_r_1, "data/processed_data/wui_r_1.rds")

Map of WUI areas (Figure 3.22; Table 3.11).

Show the code 
wui_r_1 <- terra::readRDS("data/processed_data/wui_r_1.rds")
Show the code 
f_plot <- function(.x){ plot(.x, 
                             axes = F, 
                             mar = c(0, 0, 0, 0),
                             legend = "topleft")  }

f_plot(wui_r_1$effis)
f_plot(wui_r_1$corine)
f_plot(wui_r_1$firEUrisk)
f_plot(wui_r_1$national)
(a) EFFIS
(b) CORINE
(c) FirEUrisk
(d) National databases
Figure 3.22: WUI maps of Norway based on different fuel/vegetation maps (resolution of 250 x 250 m)
Show the code 
# Summary WUI by country
wui_r_1_no <- mask(wui_r_1, no_v)
wui_r_1_se <- mask(wui_r_1, se_v)

# Tables
tbl_summary_vegetation(wui_r_1_no)
tbl_summary_vegetation(wui_r_1_se)
Table 3.11: Percentages of WUI by vegetation
(a) Norway

No WUI

WUI

Interface

Intermix

Total
corine 72.5 8.8 18.6 27.5
effis 71.4 6.9 21.8 28.6
firEUrisk 77.5 5.6 16.9 22.5
national 82.0 10.2 7.8 18.0
(b) Sweden

No WUI

WUI

Interface

Intermix

Total
corine 86.0 5.3 8.6 14.0
effis 85.8 4.9 9.3 14.2
firEUrisk 88.0 3.0 8.9 12.0
national 85.4 5.8 8.8 14.6

Plot all together

Show the code 
p_1 <- ggplot() + 
  geom_spatraster(data = wui_r) +
  facet_wrap(~lyr, nrow = 1) +
  labs(title = "d = 1.00 #/km2") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

p_3.08 <- ggplot() + 
  geom_spatraster(data = wui_r) +
  facet_wrap(~lyr, nrow = 1) +
  labs(title = "d = 3.08 #/km2") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

p_6.17 <- ggplot() + 
  geom_spatraster(data = wui_r) +
  facet_wrap(~lyr, nrow = 1) +
  labs(title = "d = 6.17 #/km2") +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5))

# NOTE: geom_spatraster change the resolution of the map
p_1 / p_3.08 / p_6.17 +
  plot_layout(guides = "collect")
Figure 3.23

4 Wildfires datasets

Show the code 
wf_no <- readxl::read_xlsx("data/wui_fires_norway_wgs84.xlsx", skip = 1) |> 
  janitor::clean_names() |> 
  st_as_sf(coords = c("geolokasjon_lengdegrad", "geolokasjon_breddegrad"),
           crs =  4326) |> 
  select(geometry) |> 
  st_transform(crs = 3035) |> 
  mutate(country = "NO")

# Infrastructures affected by wildfires in SE
wf_se_inf <- readxl::read_xlsx("data/wui_fires_sweden_epsg3006.xlsx") |> 
  janitor::clean_names() |> 
  filter(datetime >= 2016) |> 
  st_as_sf(coords = c("ecoord", "nccord"),
           crs =  3006) |> 
  select(geometry) |> 
  st_transform(crs = 3035) |> 
  mutate(country = "SE")

# Wildfires in SE (unique coordinates - one wildfire can affect several buildings)
wf_se <- distinct(wf_se_inf)

# Final dataset (NO + SE)
wf <- rbind(wf_no, wf_se)

# Map
ggplot() +
  geom_sf(data = region_sf, fill = "#FAFAFA") +
  geom_sf(data = country_sf, fill = NA) +
  geom_sf(data = wf, col = "red", size = 0.4) +
  theme_void()

We have analysed all wildfires that have affected buildings from 2016 to 2023 (i.e., 268; 88 in Norway and 180 in Sweden). The can be obtained from the Norwegian Directorate for Civil Protection (DSB) (Norway) and the Swedish Civil Contingencies Agency (MSB) (Sweden). [Note: They were not uploaded in our GitHub repository]

4.1 Wilfires and building density

Calculate the building density of the grid of 250 x 250 m where the wilfire originated.

Show the code 
ggplot() +
  tidyterra::geom_spatvector(data = region_v, fill = "white") +
  tidyterra::geom_spatraster(data = build_focal) +
  scale_fill_whitebox_c(name = "N. Buildings",
                        palette = "viridi",
                        trans = "log10") +
  geom_sf(data = wf, col = "red", size = 1) +
  theme_void()
Figure 4.1: Number of buildings per grid cells of 250 m x 250 m and wildfires (red points)
Show the code 
# Link wildfires with building density (build_focal) 
wf_build_density <- extract(build_focal, wf, touches = TRUE) |>
  # Assign 0 to NAs --> No buildings
  mutate(focal_sum = ifelse(is.na(focal_sum), 0, focal_sum)) |> 
  # Building density is the number of buildings / buffer area 
  # buffer area = circular moving window of 500 m
  mutate(bd = focal_sum / area_buffer)


ggplot(wf_build_density, aes(x = bd)) + 
  geom_histogram(bins = 50, color = "darkblue", fill = "lightblue") + 
  theme_bw()
Figure 4.2: Histogram of building density in each wildfire
Show the code 
brks <-  c(0,  1, 3, 6, 10, 20, 50, 100, 500, 1700)
wf_build_density |> 
  mutate(bd_brks = cut(bd, 
                       brks,
                       include.lowest = TRUE,
                       right = FALSE,
                       dig.lab = 4)) |> 
  group_by(bd_brks) |> 
  summarise(n = n(), 
            perc = n() / nrow(wf_build_density) * 100) |> 
  mutate(cs = cumsum(perc)) |> 
  gt() |> 
  fmt_number(columns = c("perc", "cs"),
             decimals = 1) |> 
  tab_spanner(label = md("**Wildfires**"),
              columns = c(n, perc, cs)) |> 
  cols_label(
    bd_brks = md(html("**Building**<br>**density**")),
    n  = md("**Number**"),
    perc = md("**Percentage [%]**"),
    cs = md("**Cumulative [%]**")
  ) |> 
  tab_options(
    table.font.size = 9,
    row_group.padding = px(1),
    row_group.font.weight = "bold",
    data_row.padding = px(1)
  )
Table 4.1: Number of wildfires

Building
density

Wildfires

Number

Percentage [%]

Cumulative [%]
[0,1) 62 23.1 23.1
[1,3) 13 4.9 28.0
[3,6) 3 1.1 29.1
[6,10) 11 4.1 33.2
[10,20) 31 11.6 44.8
[20,50) 50 18.7 63.4
[50,100) 40 14.9 78.4
[100,500) 39 14.6 92.9
[500,1700] 19 7.1 100.0

4.2 WUI definition

4.2.1 Effect distance to large forest

Get if the grid to which the wildfire is located was defined as WUI.

Show the code 
# Detect wui class for each point (wildfire)
extract_wui <- function(.x) { extract(.x, wf) }

wf_wui <- list(wui_effis, wui_corine, wui_national, wui_firEUrisk) |> 
  map(extract_wui) |> 
  map(as_tibble)
names(wf_wui) <- c("effis", "corine", "national", "firEUrisk")

# Tables with counts of wildfire in each WUI class
f_counts <- function(.x) {
  .x |> 
  pivot_longer(-ID) |> 
  count(name, value) |> 
  mutate(name = factor(name, 
                       levels = c("d200",
                                  "d500",
                                  "d1000",
                                  "d2000",
                                  "d2400"),
                       ordered = T)) |>
  arrange(name) |> 
  pivot_wider(names_from = name,
              values_from = n) |> 
  mutate(across(.cols = starts_with("d"),
                .fns = ~ 100 * (.x / sum(.x)),
                .names = "perc_{.col}"))
  
} 

wf_wui_counts <- wf_wui |>
  map(f_counts) |> 
  # Merge with names
  bind_rows(.id = "vegetation") 

# Print table
wf_wui_counts |> 
  # select("vegetation", "value", starts_with("perc")) |> 
  gt(groupname_col = "vegetation",
     row_group_as_column = FALSE) |> 
  fmt_number(columns = starts_with("perc"),
             decimals = 1) |> 
  cols_label(
    vegetation = "",
    value = md("**WUI class**"),
    perc_d200 = md("**d200**"),
    perc_d500 = md("**d500**"),
    perc_d1000 = md("**d1000**"),
    perc_d2000 = md("**d2000**"),
    perc_d2400 = md("**d2400**")
  )  |> 
  cols_label_with(fn = function(x) {
    janitor::make_clean_names(x, case = "title") |>
      tolower() |>
      stringr::str_replace_all("^|$", "**") |>
      md()
  }) |> 
  tab_spanner(
    label = md("**Percentage**"),
    columns = starts_with("perc")
  ) |> 
  tab_spanner(
    label = md("**Count**"),
    columns = starts_with("d")
  ) |> 
  tab_options(
    table.font.size = 9,
    row_group.background.color = "lightgrey",
    row_group.padding = px(1),
    row_group.font.weight = "bold",
    data_row.padding = px(1)
  )
Table 4.2: Percentage of wildfires in each WUI class

wui class

Count

Percentage

d200

d500

d1000

d2000

d2400

d200

d500

d1000

d2000

d2400
effis
No_wui 161 124 110 103 100 60.1 46.3 41.0 38.4 37.3
Interface 27 64 78 85 88 10.1 23.9 29.1 31.7 32.8
Intermix 80 80 80 80 80 29.9 29.9 29.9 29.9 29.9
corine
No_wui 164 130 119 108 106 61.2 48.5 44.4 40.3 39.6
Interface 35 69 80 91 93 13.1 25.7 29.9 34.0 34.7
Intermix 69 69 69 69 69 25.7 25.7 25.7 25.7 25.7
national
No_wui 180 155 148 135 129 67.2 57.8 55.2 50.4 48.1
Interface 33 58 65 78 84 12.3 21.6 24.3 29.1 31.3
Intermix 55 55 55 55 55 20.5 20.5 20.5 20.5 20.5
firEUrisk
No_wui 168 151 151 151 151 62.7 56.3 56.3 56.3 56.3
Interface 24 41 41 41 41 9.0 15.3 15.3 15.3 15.3
Intermix 76 76 76 76 76 28.4 28.4 28.4 28.4 28.4

4.2.2 Effect threshold building density

Show the code 
wui_r_1 <- terra::readRDS("data/processed_data/wui_r_1.rds")
wui_r_3 <- terra::readRDS("data/processed_data/wui_r_3.rds")
wui_r <- terra::readRDS("data/processed_data/wui_r.rds")
Show the code 
# Detect wui class for each point (wildfire)
extract_wui <- function(.x) { extract(.x, wf) }

wf_wui_bd <- list(wui_r_1, wui_r_3, wui_r) |> 
  map(extract_wui) |> 
  map(as_tibble)
names(wf_wui_bd) <- c("bd_1", "bd_3", "bd_6")

# Tables with counts og the N. wildfire in each wui class
f_counts <- function(.x) {
  .x |> 
  pivot_longer(-ID) |> 
  count(name, value) |> 
  mutate(name = factor(name, 
                       levels = c("effis",
                                  "corine",
                                  "national",
                                  "firEUrisk"),
                       ordered = T)) |>
  arrange(name) |> 
  pivot_wider(names_from = name,
              values_from = n) |> 
  mutate(across(.cols = -value,
                .fns = ~ 100 * (.x / sum(.x)),
                .names = "perc_{.col}"))
  
} 


wf_wui_bd_counts <- wf_wui_bd |>
  map(f_counts) |> 
  # Merge with names
  bind_rows(.id = "threshold") 

# Print table
wf_wui_bd_counts |> 
  mutate(definition = case_when(
    threshold == "bd_1" ~ "threshold = 1.00 #/km2",
    threshold == "bd_3" ~ "threshold = 3.08 #/km2",
    threshold == "bd_6" ~ "threshold = 6.17 #/km2"
  )) |> 
  select(-threshold) |> 
  gt(groupname_col = "definition",
     row_group_as_column = FALSE) |> 
  fmt_number(columns = starts_with("perc"),
             decimals = 1) |> 
  cols_label(
    definition = "",
    value = md("**WUI class**"),
    effis = md("**effis**"),
    corine = md("**corine**"),
    national = md("**national**"),
    firEUrisk = md("**firEUrisk**"),
    perc_effis = md("**effis**"),
    perc_corine = md("**corine**"),
    perc_national = md("**national**"),
    perc_firEUrisk = md("**firEUrisk**")
  )  |> 
  tab_spanner(
    label = md("**Percentage**"), columns = starts_with("perc")
  ) |> 
  tab_spanner(
    label = md("**Count**"), columns = c(effis, corine, national, firEUrisk)) |> 
  tab_options(
    table.font.size = 9,
    row_group.background.color = "lightgrey",
    row_group.padding = px(1),
    row_group.font.weight = "bold",
    data_row.padding = px(1)
  )
Table 4.3: Percentage of wildfires in WUI calss based on building density threshold

WUI class

Count

Percentage

effis

corine

national

firEUrisk

effis

corine

national

firEUrisk
threshold = 1.00 #/km2
No_wui 86 92 114 138 32.1 34.3 42.5 51.5
Interface 93 97 89 43 34.7 36.2 33.2 16.0
Intermix 89 79 65 87 33.2 29.5 24.3 32.5
threshold = 3.08 #/km2
No_wui 97 103 126 149 36.2 38.4 47.0 55.6
Interface 89 94 85 41 33.2 35.1 31.7 15.3
Intermix 82 71 57 78 30.6 26.5 21.3 29.1
threshold = 6.17 #/km2
No_wui 100 106 129 151 37.3 39.6 48.1 56.3
Interface 88 93 84 41 32.8 34.7 31.3 15.3
Intermix 80 69 55 76 29.9 25.7 20.5 28.4
Anderson, Hal E. 1982. “Aids to Determining Fuel Models for Estimating Fire Behavior.” INT-GTR-122. Ogden, UT: U.S. Department of Agriculture, Forest Service, Intermountain Forest; Range Experiment Station. https://doi.org/10.2737/INT-GTR-122.
Aragoneses, Elena, Mariano Garcia, and Emilio Chuvieco. 2022. “FirEUrisk_europe_fuel_map: European Fuel Map at 1 Km Resolution.” e-cienciaDatos. https://doi.org/10.21950/YABYCN.
Bar-Massada, Avi. 2021. “A Comparative Analysis of Two Major Approaches for Mapping the Wildland-Urban Interface: A Case Study in California.” Land 10 (7): 679. https://doi.org/10.3390/land10070679.
Bar-Massada, Avi, Susan I. Stewart, Roger B. Hammer, Miranda H. Mockrin, and Volker C. Radeloff. 2013. “Using Structure Locations as a Basis for Mapping the Wildland Urban Interface.” Journal of Environmental Management 128 (October): 540–47. https://doi.org/10.1016/j.jenvman.2013.06.021.
EFFIS. 2017. “European Forest Fire Information System (EFFIS) - European Fuel Map, 2017, Based on JRC Contract Number 384347 on the ‘Development of a European Fuel Map’.” European Commission.
Stewart, Susan I., Volker C. Radeloff, Roger B. Hammer, and Todd J. Hawbaker. 2007. “Defining the Wildland–Urban Interface.” Journal of Forestry 105 (4): 201–7. https://doi.org/10.1093/jof/105.4.201.