Social Vulnerability Logroño - Calculate Vulnerability

Environment

R Libraries

The relvant R libraries are imported in to the kernal:

# Load R libraries
if(!require("pacman"))
    install.packages("pacman")

p_load("sf", "tidyverse")

print("Loaded Packages:")
p_loaded()

Loading required package: pacman

[1] "Loaded Packages:"

1. 'lubridate'
2. 'forcats'
3. 'stringr'
4. 'dplyr'
5. 'purrr'
6. 'readr'
7. 'tidyr'
8. 'tibble'
9. 'ggplot2'
10. 'tidyverse'
11. 'sf'
12. 'pacman'

Output directory

# create the output directory if it does not exist
output_dir <- file.path("../..","3_outputs","Spain","Logrono","2021")
if(!dir.exists(output_dir)){
    dir.create(output_dir, recursive = TRUE)
    print(paste0(output_dir, " created"))
}

Set the GUID

GUID <- c("CCA", "CPRO", "CMUN", "CDIS", "CSEC")
GUID_census_indicator_data <- c("ccaa", "CPRO", "CMUN", "dist", "secc")
GUID_length = 5

Load Data

Import the data

# Load census data
census_indicator_data <- read.csv("../../2_pipeline/Spain/Logrono/1a_CensusData/2021/censusDataZ.csv")
# Update the census indicator data ID to the common GUID
colnames(census_indicator_data)[colnames(census_indicator_data) %in% GUID_census_indicator_data] = GUID
# Change all the GUID column data type to integer, this helps to merge datasets later in the code
census_indicator_data[GUID] <- lapply(census_indicator_data[GUID], as.integer)
# Print head of census indicator data
head(census_indicator_data)

# Load Coperncius data: tree cover density (TCD) and imperviousness density (IMP)
# Both these datasets already uses the GUID
# Both these datasets contain the census area geometry, so we extract from one of the datasets for use later in the code
tcd_indicator_data <- st_read("../../2_pipeline/Spain/1b_Copernicus/2021/census_areas_TCD.geojson")
imd_indicator_data <- st_read("../../2_pipeline/Spain/1b_Copernicus/2021/census_areas_IMD.geojson")
# Get the geospatial data from the TCD data (the IMD data also has same spatial data)
oa <- subset(tcd_indicator_data, select = c(GUID, 'geometry'))
# Remove the geospatial geometry from the TCS and IMD datasets
tcd_indicator_data <- st_drop_geometry(tcd_indicator_data)
imd_indicator_data <- st_drop_geometry(imd_indicator_data)
# Change all the GUID column data type to integer, this helps to merge datasets later in the code
imd_indicator_data[GUID] <- lapply(imd_indicator_data[GUID], as.integer)
tcd_indicator_data[GUID] <- lapply(tcd_indicator_data[GUID], as.integer)
oa[GUID] <- lapply(st_drop_geometry(oa[GUID]), as.integer)

# Load vulnerability mapping information from the config file
## This mapping information is used to help guide the amalgamation of the data.
## Weighting can be changed in this file, depending on the scenario.
## Scenario 1 (best case scenario): Weighting values 1 or -1:
##  where 1 means no change
##  or -1 means all the indicator values are multiplied by -1, resulting in an inverse indicator.
## Scenario 2: Weighting values 0.5 or -0.5:
##  For domains with just a single indicators or where there is a lack of information related to missing indicators. 
##  For these domains the weights are halved using a weight of 0.5, or -0.5 for an inverse indicator.
##  Therefore the influence of these indicators are reduced in half.
## Other scenarios are supported by using other decimal numbers if decided for a particular dataset.
indicator_mapping <- read.csv("config/vulnerabilityIndicatorMappings.csv", header=TRUE, sep=",", stringsAsFactors = FALSE, fileEncoding="UTF-8-BOM")

# Print up to 100 rows of vulnerabiltiy mapping config file
head(indicator_mapping,100)

A data.frame: 6 × 20

	CCA	CPRO	CMUN	CDIS	CSEC	early_childhood_boy	early_childhood_girl	age_middle_to_oldest_old_male	age_middle_to_oldest_old_female	disability	one_parent_households	dependants	unemployment	attending_university	no_higher_education	foreign_nationals	rented	primary_school_age	one_person_households	year_built
	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	17	26	1	1	1	-0.19258164	-0.1467916	0.4818082	-0.10021618	1.2753218	-1.459709836	-0.0938240	-1.00214549	1.2303104	0.07143877	0.1446865	-0.5073936	0.30936208	1.0319635	-0.94812767
2	17	26	2	1	1	0.23745377	0.3146154	-0.2611307	-0.31519260	0.3880289	-0.007254482	0.1987363	0.87022439	0.2231589	0.27339067	0.7793089	-0.5005024	0.16550926	0.4903007	0.24638629
3	17	26	3	1	1	-0.96122760	-0.1140934	0.3255496	0.83585043	-0.4992639	-1.459709836	-0.0938240	-0.37802220	0.7267347	0.20607337	0.3985355	-1.0793008	0.26195604	0.7214991	2.13797352
4	17	26	4	1	1	0.01083739	-1.4002199	0.5976241	-1.28188374	-1.3865567	NA	-1.5566254	-1.62626879	-1.2875683	1.61973671	-1.2514828	NA	-1.57053273	NA	NA
5	17	26	5	1	1	-0.10396496	0.1640514	-0.1893316	0.06098884	1.2753218	-0.118741371	0.3450164	1.18228604	-0.2804168	0.27339067	0.6523844	-0.6041755	0.05555571	-0.1104731	-0.08916627
6	17	26	5	2	1	2.13595981	1.4136596	-0.7006675	-0.72507994	0.3880289	1.157623344	1.5152575	-0.06596055	0.2231589	-0.06319583	0.1446865	-0.3682738	1.19901504	-0.4718037	-0.42805727

Reading layer `census_areas_TCD' from data source 
  `/Cities/2_pipeline/Spain/1b_Copernicus/2021/census_areas_TCD.geojson' 
  using driver `GeoJSON'
Simple feature collection with 36333 features and 22 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -1004502 ymin: 3132130 xmax: 1126932 ymax: 4859240
Projected CRS: ETRS89 / UTM zone 30N
Reading layer `census_areas_IMD' from data source 
  `/Cities/2_pipeline/Spain/1b_Copernicus/2021/census_areas_IMD.geojson' 
  using driver `GeoJSON'
Simple feature collection with 36333 features and 22 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -1004502 ymin: 3132130 xmax: 1126932 ymax: 4859240
Projected CRS: ETRS89 / UTM zone 30N

A data.frame: 17 × 9

	domain	indicator	sensitivity	prepare	respond	recover	adaptive_capacity	enhanced_exposure	weight
	<chr>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>
1	age	early_childhood_boy	1	0	0	0	0	0	1
2	age	early_childhood_girl	1	0	0	0	0	0	1
3	age	age_middle_to_oldest_old_male	1	0	0	0	0	0	1
4	age	age_middle_to_oldest_old_female	1	0	0	0	0	0	1
5	health	disability	1	0	0	0	0	0	1
6	income	one_parent_households	0	1	1	1	1	0	1
7	income	dependants	0	1	1	1	1	0	1
8	income	unemployment	0	1	1	1	1	0	1
9	income	attending_university	0	1	1	1	1	0	1
10	info_access_use	no_higher_education	0	1	1	1	1	0	1
11	local_knowledge	foreign_nationals	0	1	1	0	1	0	1
12	tenure	rented	0	1	0	0	1	0	1
13	social_network	primary_school_age	0	0	1	1	1	0	-1
14	social_network	one_person_households	0	0	1	1	1	0	1
15	housing_characteristics	year_built	0	0	0	0	0	1	1
16	physical_environment	impervious	0	0	0	0	0	1	1
17	physical_environment	tree_cover_density	0	0	0	0	0	1	1

Prepare Data

Combine data into a single indicator dataset

# combine census data with copernicus TCD and IMD data (without geospatial data to avoid duplication)
indicator_data <- merge(census_indicator_data, tcd_indicator_data, by=GUID, all.x=TRUE)
indicator_data <- merge(indicator_data, imd_indicator_data, by=GUID, all.x=TRUE)

# trim the columns
indicator_data <- subset(indicator_data, select=c(names(census_indicator_data), 'tree_cover_density', 'impervious'))

# Set missing data fields to zero (0)
indicator_data[is.na(indicator_data)] <- 0

# Print the first part of the indicators, which are now collated into one table
head(indicator_data)

A data.frame: 6 × 22

	CCA	CPRO	CMUN	CDIS	CSEC	early_childhood_boy	early_childhood_girl	age_middle_to_oldest_old_male	age_middle_to_oldest_old_female	disability	⋯	unemployment	attending_university	no_higher_education	foreign_nationals	rented	primary_school_age	one_person_households	year_built	tree_cover_density	impervious
	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	17	26	1	1	1	-0.19258164	-0.1467916	0.48180824	-0.10021618	1.2753218	⋯	-1.0021455	1.2303104	0.071438769	0.1446865	-0.5073936	0.3093621	1.0319635	-0.9481277	-0.4793341	-1.232806
2	17	26	10	1	1	-0.36188692	-0.3284478	1.32081073	0.44848018	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	-0.6060650	0.0000000	0.0000000	0.6403542	-1.105968
3	17	26	100	1	1	0.01083739	-1.4002199	0.23798535	-0.52201032	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	-0.9083608	0.0000000	0.0000000	0.5002384	-1.264703
4	17	26	101	1	1	-1.36078806	-1.4002199	-0.22155305	0.07955613	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	-1.5705327	0.0000000	0.0000000	-1.7574746	-1.265740
5	17	26	102	1	1	0.10257030	-0.4581527	0.05129390	0.59299501	0.3880289	⋯	1.0262552	0.2231589	0.138756071	-0.1091625	-0.7479582	-0.4402011	0.3954425	1.2637707	0.4972984	-1.255886
6	17	26	102	1	2	0.64987009	0.6302157	-0.04065368	0.64469364	0.3880289	⋯	0.7141936	-0.7839926	0.004121467	0.1446865	-0.9688721	0.7134064	0.8394153	0.7211874	0.4394261	-1.178694

Weight the indicator datas

# Get the indicator weighting, previously loaded from the config file
indicator_weighting <- indicator_mapping %>% select('indicator', 'weight')
indicator_weighting <- indicator_weighting %>% spread(key = 'indicator', value = 'weight')

# Get the column names and weights
names <- names(indicator_weighting)
weights <- indicator_weighting[, names]

# Copy and rename the dataset
indicator_data_weighted <- indicator_data
head(indicator_data_weighted) 

# Multiply the indicators by the config file weighting
indicator_data_weighted[, names] <- sweep(indicator_data_weighted[, names], 2, unlist(weights[, names]), "*")
head(indicator_data_weighted)

A data.frame: 6 × 22

	CCA	CPRO	CMUN	CDIS	CSEC	early_childhood_boy	early_childhood_girl	age_middle_to_oldest_old_male	age_middle_to_oldest_old_female	disability	⋯	unemployment	attending_university	no_higher_education	foreign_nationals	rented	primary_school_age	one_person_households	year_built	tree_cover_density	impervious
	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	17	26	1	1	1	-0.19258164	-0.1467916	0.48180824	-0.10021618	1.2753218	⋯	-1.0021455	1.2303104	0.071438769	0.1446865	-0.5073936	0.3093621	1.0319635	-0.9481277	-0.4793341	-1.232806
2	17	26	10	1	1	-0.36188692	-0.3284478	1.32081073	0.44848018	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	-0.6060650	0.0000000	0.0000000	0.6403542	-1.105968
3	17	26	100	1	1	0.01083739	-1.4002199	0.23798535	-0.52201032	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	-0.9083608	0.0000000	0.0000000	0.5002384	-1.264703
4	17	26	101	1	1	-1.36078806	-1.4002199	-0.22155305	0.07955613	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	-1.5705327	0.0000000	0.0000000	-1.7574746	-1.265740
5	17	26	102	1	1	0.10257030	-0.4581527	0.05129390	0.59299501	0.3880289	⋯	1.0262552	0.2231589	0.138756071	-0.1091625	-0.7479582	-0.4402011	0.3954425	1.2637707	0.4972984	-1.255886
6	17	26	102	1	2	0.64987009	0.6302157	-0.04065368	0.64469364	0.3880289	⋯	0.7141936	-0.7839926	0.004121467	0.1446865	-0.9688721	0.7134064	0.8394153	0.7211874	0.4394261	-1.178694

A data.frame: 6 × 22

	CCA	CPRO	CMUN	CDIS	CSEC	early_childhood_boy	early_childhood_girl	age_middle_to_oldest_old_male	age_middle_to_oldest_old_female	disability	⋯	unemployment	attending_university	no_higher_education	foreign_nationals	rented	primary_school_age	one_person_households	year_built	tree_cover_density	impervious
	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	17	26	1	1	1	-0.19258164	-0.1467916	0.48180824	-0.10021618	1.2753218	⋯	-1.0021455	1.2303104	0.071438769	0.1446865	-0.5073936	-0.3093621	1.0319635	-0.9481277	-0.4793341	-1.232806
2	17	26	10	1	1	-0.36188692	-0.3284478	1.32081073	0.44848018	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	0.6060650	0.0000000	0.0000000	0.6403542	-1.105968
3	17	26	100	1	1	0.01083739	-1.4002199	0.23798535	-0.52201032	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	0.9083608	0.0000000	0.0000000	0.5002384	-1.264703
4	17	26	101	1	1	-1.36078806	-1.4002199	-0.22155305	0.07955613	-1.3865567	⋯	-1.6262688	-1.2875683	1.619736710	-1.2514828	0.0000000	1.5705327	0.0000000	0.0000000	-1.7574746	-1.265740
5	17	26	102	1	1	0.10257030	-0.4581527	0.05129390	0.59299501	0.3880289	⋯	1.0262552	0.2231589	0.138756071	-0.1091625	-0.7479582	0.4402011	0.3954425	1.2637707	0.4972984	-1.255886
6	17	26	102	1	2	0.64987009	0.6302157	-0.04065368	0.64469364	0.3880289	⋯	0.7141936	-0.7839926	0.004121467	0.1446865	-0.9688721	-0.7134064	0.8394153	0.7211874	0.4394261	-1.178694

# Histogram visualisation of weighted indicators
index_start = GUID_length+1
index_end = ncol(indicator_data_weighted)
indicator_columns <- colnames(indicator_data_weighted)[index_start:index_end]
for( current_indicator_column in indicator_columns ) {
    indicator_filtered <- indicator_data_weighted[,current_indicator_column] 
    indicator_filtered[indicator_filtered == "NaN"] <- 0

    title <- paste("'", current_indicator_column, "' domain histogram", sep = "")
    x_label <- paste("'", current_indicator_column, "' z-score", sep = "")
    y_label <- paste("Count", sep = "")
    hist(indicator_filtered, breaks="FD", col="grey", labels = FALSE, main=title, xlab=x_label, ylab=y_label)
    box("figure", lwd = 4)
}

Process social vulnerability scores

Calculate domain scores

# Get the domains and their associated indicator ID
domain_indicators <- indicator_mapping %>% select('domain', 'indicator')

# Get a vector/array of the unique domain names
unique_domains <- unique(domain_indicators$domain)

# Initialise the domain score dataset with the GUID
domain_scores <- indicator_data_weighted %>% select(all_of(GUID))

# Loop through each domain
for (current_domain in unique_domains) {
    # Identify which indicators are used within this domain (current_domain)
    current_domain_info <- domain_indicators %>% filter(domain == current_domain)

    # Count the number of indicators in this domain
    domain_indicator_count <- length(current_domain_info$indicator)

    # Get a vector/array of the indicators used by this domain, and add the GUID column name
    current_domain_indicators <- current_domain_info$indicator
    current_domain_indicators <- (c(GUID, current_domain_indicators))

    # filter the dataset to only use the indicators in the domain
    current_domain_data <- indicator_data_weighted[current_domain_indicators]

    # Calculate the internal weight distribution for the indicators within this domain,
    # using an equal weight distribution across this domain
    internal_domain_weight <- 1.0 / domain_indicator_count

    # Internally weight the data for this domain
    current_domain_data_weighted <- current_domain_data %>% mutate_if(is.double, function(x) {x*internal_domain_weight})

    # Sum each data row to get the total score for the domain
    index_start = GUID_length+1
    index_end = domain_indicator_count+5
    current_domain_data_weighted[, current_domain] <- rowSums(current_domain_data_weighted[index_start:index_end], na.rm = TRUE)

    # Add the current domain score to the overall results
    domain_indicator_score <- current_domain_data_weighted %>% select(all_of(GUID), all_of(current_domain))
    domain_scores <- merge(domain_scores, domain_indicator_score, by=GUID)
}

# Print the first part of the domain z-scores, which are now collated into one table
head(domain_scores)

A data.frame: 6 × 14

	CCA	CPRO	CMUN	CDIS	CSEC	age	health	income	info_access_use	local_knowledge	tenure	social_network	housing_characteristics	physical_environment
	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	17	26	1	1	1	0.01055472	1.2753218	-0.3313422	0.071438769	0.1446865	-0.5073936	0.36130070	-0.9481277	-0.8560700
2	17	26	10	1	1	0.26973904	-1.3865567	-1.1176156	1.619736710	-1.2514828	0.0000000	0.30303248	0.0000000	-0.2328068
3	17	26	100	1	1	-0.41835186	-1.3865567	-1.1176156	1.619736710	-1.2514828	0.0000000	0.45418041	0.0000000	-0.3822322
4	17	26	101	1	1	-0.72575121	-1.3865567	-1.1176156	1.619736710	-1.2514828	0.0000000	0.78526637	0.0000000	-1.5116073
5	17	26	102	1	1	0.07217662	0.3880289	0.1439333	0.138756071	-0.1091625	-0.7479582	0.41782180	1.2637707	-0.3792937
6	17	26	102	1	2	0.47103144	0.3880289	0.2297142	0.004121467	0.1446865	-0.9688721	0.06300444	0.7211874	-0.3696338

# Histogram visualisation of domain z-scores
for( current_domain in unique_domains ) {
    domain_scores_filtered <- domain_scores[,current_domain] 
    domain_scores_filtered[domain_scores_filtered == "NaN"] <- 0

    title <- paste("'", current_domain, "' domain histogram", sep = "")
    x_label <- paste("'", current_domain, "' z-score", sep = "")
    y_label <- paste("Count", sep = "")
    hist(domain_scores_filtered, breaks="FD", col="grey", labels=FALSE, main=title, xlab=x_label, ylab=y_label)
    box("figure", lwd = 4)
}

Calculate dimension scores

Need to collate the domains into the dimensions

# Create a vector/array of the dimension names
dimensions <- c('sensitivity', 'prepare', 'respond', 'recover', 'adaptive_capacity', 'enhanced_exposure')

# Get the dimension and their associated indicator ID
dimension_indicators <- indicator_mapping %>% select(c('domain', all_of(dimensions)))
head(dimension_indicators)

# Initialise the dimensions score dataset with the GUID
dimension_scores <- indicator_data_weighted %>% select(all_of(GUID))

# loop through each of the dimensions and:
for (current_dimension in dimensions){
    # Identify which indicators are used within this dimension (current_dimension)
    # Then select the indicators marked with value 1, which means that the indicator is part of the dimension
    current_dimension_info <- dimension_indicators %>% select(c('domain', all_of(current_dimension)))
    current_dimension_info <- current_dimension_info %>% filter(dimension_indicators[, current_dimension] == 1)

    # Get a array/vector of the unique domains in this dimension
    current_dimension_domains <- unique(current_dimension_info$domain)

    # Count the number of domains in this dimension
    dimension_domain_count <- length(current_dimension_domains)

    # Filter the domain scores dataset to only use the domains in the dimension, and add the GUID column name
    current_dimension_data <- domain_scores %>% select(c(all_of(GUID), all_of(current_dimension_domains)))  

    # Sum each data row to get the total score for the dimension
    index_start = GUID_length+1
    index_end = dimension_domain_count+5
    current_dimension_data[, current_dimension] <- rowSums(current_dimension_data[index_start:index_end], na.rm = TRUE)

    # Add the current dimension score to the overall results
    dimension_indicator_score <- current_dimension_data %>% select(all_of(GUID), all_of(current_dimension))
    dimension_scores <- merge(dimension_scores, dimension_indicator_score, by=GUID)  
}

# generate z-scores with the scale function in order to standardise the dimension data
dimension_scores <- dimension_scores %>% mutate_if(is.double, scale)

# Print the first part of the dimension scores, which are now collated into one table
head(dimension_scores)

A data.frame: 6 × 7

	domain	sensitivity	prepare	respond	recover	adaptive_capacity	enhanced_exposure
	<chr>	<int>	<int>	<int>	<int>	<int>	<int>
1	age	1	0	0	0	0	0
2	age	1	0	0	0	0	0
3	age	1	0	0	0	0	0
4	age	1	0	0	0	0	0
5	health	1	0	0	0	0	0
6	income	0	1	1	1	1	0

A data.frame: 6 × 11

	CCA	CPRO	CMUN	CDIS	CSEC	sensitivity	prepare	respond	recover	adaptive_capacity	enhanced_exposure
	<int>	<int>	<int>	<int>	<int>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>
1	17	26	1	1	1	1.1377171	-0.3378133	0.18457796	0.1073270	-0.15329655	-1.00827077
2	17	26	10	1	1	-0.9881374	-0.4065854	-0.33477437	0.8522391	-0.26183757	0.04893220
3	17	26	100	1	1	-1.5969460	-0.4065854	-0.22140412	1.0122262	-0.17316713	-0.05159845
4	17	26	101	1	1	-1.8689267	-0.4065854	0.02693072	1.3626742	0.02106336	-0.81142134
5	17	26	102	1	1	0.4071804	-0.3116725	0.44354790	0.7414772	-0.09187444	0.80062025
6	17	26	102	1	2	0.7600790	-0.3203095	0.33117214	0.3141994	-0.30936550	0.44207919

# Histogram visualisation of dimension z-scores
for( current_dimension in dimensions ){
    dimension_scores_filtered <- dimension_scores[,current_dimension] 
    dimension_scores_filtered[dimension_scores_filtered == "NaN"] <- 0
   
    title <- paste("'", current_dimension, "' dimension histogram", sep = "")
    x_label <- paste("'", current_dimension, "' z-score", sep = "")
    y_label <- paste("Count", sep = "")
    hist(dimension_scores_filtered, breaks="FD", col="grey", labels=FALSE, main=title, xlab=x_label, ylab=y_label)
    box("figure", lwd = 4)
}

Calculate vulnerability score

# Initialise the vulnerability score dataset with the GUID
vulnerability_scores <- domain_scores %>% select(all_of(GUID))

#sum the domains to create a total overall score of vulnerability
index_start = GUID_length+1
index_end = ncol(domain_scores)
vulnerability_scores$social_vulnerability <- rowSums(domain_scores[index_start:index_end], na.rm = TRUE)

# generate z-scores with the scale function in order to standardise the vulnerability data
vulnerability_scores <- vulnerability_scores %>% mutate_if(is.double, scale)

# Print the first part of the vulnerability scores, which are now collated into one table
head(vulnerability_scores)

A data.frame: 6 × 6

	CCA	CPRO	CMUN	CDIS	CSEC	social_vulnerability
	<int>	<int>	<int>	<int>	<int>	<dbl[,1]>
1	17	26	1	1	1	-0.1413894
2	17	26	10	1	1	-0.4444888
3	17	26	100	1	1	-0.6491854
4	17	26	101	1	1	-0.9789365
5	17	26	102	1	1	0.4454417
6	17	26	102	1	2	0.2948932

# Histogram visualisation of Vulnerability z-scores
title <- paste("'Vulnerability' histogram", sep = "")
x_label <- paste("'Vulnerability' z-score", sep = "")
y_label <- paste("Count", sep = "")
hist(vulnerability_scores$social_vulnerability, breaks="FD", col="grey", labels=FALSE, main=title, xlab=x_label, ylab=y_label)
box("figure", lwd = 4)

# Merge all the indicators, domains, dimensions, and total vulnerability into one dataset
output_dataset <- merge(indicator_data_weighted, domain_scores, by=GUID)
output_dataset <- merge(output_dataset, dimension_scores, by=GUID)
output_dataset <- merge(output_dataset, vulnerability_scores, by=GUID)

head(output_dataset)

A data.frame: 6 × 38

	CCA	CPRO	CMUN	CDIS	CSEC	early_childhood_boy	early_childhood_girl	age_middle_to_oldest_old_male	age_middle_to_oldest_old_female	disability	⋯	social_network	housing_characteristics	physical_environment	sensitivity	prepare	respond	recover	adaptive_capacity	enhanced_exposure	social_vulnerability
	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>
1	17	26	1	1	1	-0.19258164	-0.1467916	0.48180824	-0.10021618	1.2753218	⋯	0.36130070	-0.9481277	-0.8560700	1.1377171	-0.3378133	0.18457796	0.1073270	-0.15329655	-1.00827077	-0.1413894
2	17	26	10	1	1	-0.36188692	-0.3284478	1.32081073	0.44848018	-1.3865567	⋯	0.30303248	0.0000000	-0.2328068	-0.9881374	-0.4065854	-0.33477437	0.8522391	-0.26183757	0.04893220	-0.4444888
3	17	26	100	1	1	0.01083739	-1.4002199	0.23798535	-0.52201032	-1.3865567	⋯	0.45418041	0.0000000	-0.3822322	-1.5969460	-0.4065854	-0.22140412	1.0122262	-0.17316713	-0.05159845	-0.6491854
4	17	26	101	1	1	-1.36078806	-1.4002199	-0.22155305	0.07955613	-1.3865567	⋯	0.78526637	0.0000000	-1.5116073	-1.8689267	-0.4065854	0.02693072	1.3626742	0.02106336	-0.81142134	-0.9789365
5	17	26	102	1	1	0.10257030	-0.4581527	0.05129390	0.59299501	0.3880289	⋯	0.41782180	1.2637707	-0.3792937	0.4071804	-0.3116725	0.44354790	0.7414772	-0.09187444	0.80062025	0.4454417
6	17	26	102	1	2	0.64987009	0.6302157	-0.04065368	0.64469364	0.3880289	⋯	0.06300444	0.7211874	-0.3696338	0.7600790	-0.3203095	0.33117214	0.3141994	-0.30936550	0.44207919	0.2948932

Correlations

# check the correlations
correlation <- cor(output_dataset %>% select(-c(all_of(GUID))), use="pairwise.complete.obs")
correlation

A matrix: 33 × 33 of type dbl

	early_childhood_boy	early_childhood_girl	age_middle_to_oldest_old_male	age_middle_to_oldest_old_female	disability	one_parent_households	dependants	unemployment	attending_university	no_higher_education	⋯	social_network	housing_characteristics	physical_environment	sensitivity	prepare	respond	recover	adaptive_capacity	enhanced_exposure	social_vulnerability
early_childhood_boy	1.00000000	0.759540053	-0.44979947	-0.372631631	0.23018526	0.33038590	0.67057589	0.35614054	0.23694332	-0.36085632	⋯	-0.548326585	0.1050579664	0.35289264	0.41094687	0.34499285	0.020408818	-0.42216457	0.12155807	0.29297386	0.330183124
early_childhood_girl	0.75954005	1.000000000	-0.46707295	-0.365530974	0.26364051	0.33986362	0.69435577	0.36004951	0.29591717	-0.39610178	⋯	-0.596400846	0.0911557556	0.38809959	0.43829722	0.35423522	0.003914884	-0.47404722	0.10950479	0.30735942	0.339651432
age_middle_to_oldest_old_male	-0.44979947	-0.467072951	1.00000000	0.554816775	-0.29036166	-0.45598571	-0.62497380	-0.34357715	-0.42706479	0.43854968	⋯	0.560920526	-0.0976966424	-0.36515398	-0.11579601	-0.32214094	0.065933061	0.44983077	-0.09107455	-0.29644999	-0.216741127
age_middle_to_oldest_old_female	-0.37263163	-0.365530974	0.55481678	1.000000000	-0.06347671	-0.45546676	-0.45061486	-0.01977655	-0.24756632	0.18899810	⋯	0.579507014	0.0464677389	-0.01493317	0.12447682	-0.16591374	0.204713854	0.38051051	0.08636588	0.01859601	0.094105996
disability	0.23018526	0.263640512	-0.29036166	-0.063476710	1.00000000	-0.04926880	0.44524799	0.66041158	0.24108336	-0.35998489	⋯	-0.104345540	0.1140625145	0.31311428	0.91574394	0.24743845	0.269148367	-0.12279689	0.21968497	0.27264011	0.541205238
one_parent_households	0.33038590	0.339863620	-0.45598571	-0.455466763	-0.04926880	1.00000000	0.44098288	0.08855721	0.41209063	-0.20215971	⋯	-0.633596005	0.2369424946	0.42875198	-0.09694488	0.50291888	-0.072330693	-0.24962760	0.25320804	0.42248985	0.283327135
dependants	0.67057589	0.694355771	-0.62497380	-0.450614861	0.44524799	0.44098288	1.00000000	0.56720146	0.55171806	-0.66446988	⋯	-0.756633822	0.0829294936	0.51345917	0.45794746	0.38170650	-0.050928458	-0.66309316	0.06572621	0.38370731	0.357862941
unemployment	0.35614054	0.360049506	-0.34357715	-0.019776546	0.66041158	0.08855721	0.56720146	1.00000000	0.33617823	-0.41956319	⋯	-0.201583893	0.2543291753	0.52983989	0.66236418	0.51846520	0.465082702	-0.08601945	0.46813351	0.49867787	0.682300003
attending_university	0.23694332	0.295917170	-0.42706479	-0.247566321	0.24108336	0.41209063	0.55171806	0.33617823	1.00000000	-0.80231522	⋯	-0.444542590	-0.0998329821	0.45902384	0.18194667	0.08596680	-0.311683928	-0.61881770	-0.11091324	0.23712302	0.110056131
no_higher_education	-0.36085632	-0.396101779	0.43854968	0.188998102	-0.35998489	-0.20215971	-0.66446988	-0.41956319	-0.80231522	1.00000000	⋯	0.359087857	0.1816427966	-0.46586422	-0.34713213	0.08891632	0.439654954	0.80803001	0.26081266	-0.19176042	-0.069422773
foreign_nationals	0.42604982	0.453075887	-0.33707442	-0.086558180	0.47484750	0.13940285	0.55855872	0.70132633	0.16908401	-0.17722844	⋯	-0.344220134	0.3930076266	0.47504867	0.52088805	0.76986598	0.706907540	-0.06089854	0.67454294	0.54754029	0.761204213
rented	0.19893150	0.200527443	-0.26872461	-0.138932327	0.01728583	0.58359421	0.19886232	0.19663693	0.25030330	-0.15646775	⋯	-0.249652474	0.2445553823	0.54975446	0.01348079	0.68961198	0.128269631	-0.06156210	0.63113802	0.50565219	0.549539325
primary_school_age	-0.65679028	-0.690484814	0.55898237	0.504557789	-0.28672716	-0.40852694	-0.84126926	-0.35945994	-0.38544019	0.50117885	⋯	0.839691152	-0.0764614131	-0.39566986	-0.31645101	-0.31384459	0.191164019	0.70772433	0.04573707	-0.30332716	-0.217873743
one_person_households	-0.22156559	-0.267393210	0.35142193	0.443686124	0.13659002	-0.64327666	-0.37761112	0.04895747	-0.34213141	0.06656413	⋯	0.799861106	-0.0008684821	-0.25593992	0.18857068	-0.31411959	0.180859976	0.40075833	0.02717408	-0.16662814	0.003512541
year_built	0.10505797	0.091155756	-0.09769664	0.046467739	0.11406251	0.23694249	0.08292949	0.25432918	-0.09983298	0.18164280	⋯	-0.049408896	1.0000000000	0.26356151	0.13299007	0.49313913	0.486775885	0.27094304	0.51053624	0.77979729	0.650035629
tree_cover_density	0.38502113	0.414558106	-0.31956017	-0.025669031	0.29130747	0.29541420	0.53298471	0.42060960	0.34500674	-0.38288149	⋯	-0.428694837	0.1803864548	0.87274526	0.35824276	0.36450898	0.087039689	-0.34561834	0.19752132	0.67620326	0.520913738
impervious	0.20875349	0.239541111	-0.30448896	0.001689509	0.24169082	0.44708573	0.33428776	0.49015701	0.44640218	-0.41631699	⋯	-0.247471063	0.2761895237	0.83655012	0.24602585	0.55573889	0.131326510	-0.20208813	0.48739139	0.71103367	0.645888772
age	0.51441045	0.508826242	0.35018899	0.448288787	0.07684377	-0.13240437	0.15883005	0.19368368	-0.07782264	-0.07103765	⋯	-0.002360354	0.0795870159	0.19811287	0.47094367	0.11592014	0.161919237	-0.03615852	0.12425339	0.17701967	0.300376066
health	0.23018526	0.263640512	-0.29036166	-0.063476710	1.00000000	-0.04926880	0.44524799	0.66041158	0.24108336	-0.35998489	⋯	-0.104345540	0.1140625145	0.31311428	0.91574394	0.24743845	0.269148367	-0.12279689	0.21968497	0.27264011	0.541205238
income	0.53809505	0.570891186	-0.62267021	-0.389151337	0.44841072	0.63584339	0.86707123	0.68303775	0.77851407	-0.71256967	⋯	-0.680471207	0.1555872409	0.65093009	0.41823692	0.49628699	0.012749347	-0.54883841	0.22454098	0.51715347	0.484368003
info_access_use	-0.36085632	-0.396101779	0.43854968	0.188998102	-0.35998489	-0.20215971	-0.66446988	-0.41956319	-0.80231522	1.00000000	⋯	0.359087857	0.1816427966	-0.46586422	-0.34713213	0.08891632	0.439654954	0.80803001	0.26081266	-0.19176042	-0.069422773
local_knowledge	0.42604982	0.453075887	-0.33707442	-0.086558180	0.47484750	0.13940285	0.55855872	0.70132633	0.16908401	-0.17722844	⋯	-0.344220134	0.3930076266	0.47504867	0.52088805	0.76986598	0.706907540	-0.06089854	0.67454294	0.54754029	0.761204213
tenure	0.19893150	0.200527443	-0.26872461	-0.138932327	0.01728583	0.58359421	0.19886232	0.19663693	0.25030330	-0.15646775	⋯	-0.249652474	0.2445553823	0.54975446	0.01348079	0.68961198	0.128269631	-0.06156210	0.63113802	0.50565219	0.549539325
social_network	-0.54832658	-0.596400846	0.56092053	0.579507014	-0.10434554	-0.63359601	-0.75663382	-0.20158389	-0.44454259	0.35908786	⋯	1.000000000	-0.0494088961	-0.40120270	-0.09327390	-0.38253686	0.226944393	0.68461214	0.04498125	-0.29044958	-0.137323581
housing_characteristics	0.10505797	0.091155756	-0.09769664	0.046467739	0.11406251	0.23694249	0.08292949	0.25432918	-0.09983298	0.18164280	⋯	-0.049408896	1.0000000000	0.26356151	0.13299007	0.49313913	0.486775885	0.27094304	0.51053624	0.77979729	0.650035629
physical_environment	0.35289264	0.388099588	-0.36515398	-0.014933172	0.31311428	0.42875198	0.51345917	0.52983989	0.45902384	-0.46586422	⋯	-0.401202696	0.2635615089	1.00000000	0.35686743	0.53128199	0.126104958	-0.32487143	0.39044819	0.80942177	0.677580930
sensitivity	0.41094687	0.438297216	-0.11579601	0.124476818	0.91574394	-0.09694488	0.45794746	0.66236418	0.18194667	-0.34713213	⋯	-0.093273903	0.1329900657	0.35686743	1.00000000	0.26563892	0.303382929	-0.12321838	0.24444113	0.31255713	0.599884974
prepare	0.34499285	0.354235222	-0.32214094	-0.165913736	0.24743845	0.50291888	0.38170650	0.51846520	0.08596680	0.08891632	⋯	-0.382536861	0.4931391345	0.53128199	0.26563892	1.00000000	0.690091341	0.15896290	0.90579806	0.64498971	0.835927532
respond	0.02040882	0.003914884	0.06593306	0.204713854	0.26914837	-0.07233069	-0.05092846	0.46508270	-0.31168393	0.43965495	⋯	0.226944393	0.4867758850	0.12610496	0.30338293	0.69009134	1.000000000	0.66294354	0.85021883	0.37816498	0.702116874
recover	-0.42216457	-0.474047225	0.44983077	0.380510515	-0.12279689	-0.24962760	-0.66309316	-0.08601945	-0.61881770	0.80803001	⋯	0.684612142	0.2709430415	-0.32487143	-0.12321838	0.15896290	0.662943544	1.00000000	0.48583097	-0.04589749	0.185101540
adaptive_capacity	0.12155807	0.109504790	-0.09107455	0.086365882	0.21968497	0.25320804	0.06572621	0.46813351	-0.11091324	0.26081266	⋯	0.044981247	0.5105362367	0.39044819	0.24444113	0.90579806	0.850218830	0.48583097	1.00000000	0.56418215	0.840851042
enhanced_exposure	0.29297386	0.307359416	-0.29644999	0.018596006	0.27264011	0.42248985	0.38370731	0.49867787	0.23712302	-0.19176042	⋯	-0.290449581	0.7797972928	0.80942177	0.31255713	0.64498971	0.378164977	-0.04589749	0.56418215	1.00000000	0.835445536
social_vulnerability	0.33018312	0.339651432	-0.21674113	0.094105996	0.54120524	0.28332714	0.35786294	0.68230000	0.11005613	-0.06942277	⋯	-0.137323581	0.6500356288	0.67758093	0.59988497	0.83592753	0.702116874	0.18510154	0.84085104	0.83544554	1.000000000

Add geometry

# add st_drop_geometry
output_dataset_geom <- merge(output_dataset, oa, by.x=GUID, by.y=GUID, all.x = TRUE)
head(output_dataset_geom)

A data.frame: 6 × 39

	CCA	CPRO	CMUN	CDIS	CSEC	early_childhood_boy	early_childhood_girl	age_middle_to_oldest_old_male	age_middle_to_oldest_old_female	disability	⋯	housing_characteristics	physical_environment	sensitivity	prepare	respond	recover	adaptive_capacity	enhanced_exposure	social_vulnerability	geometry
	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<dbl[,1]>	<MULTIPOLYGON [m]>
1	17	26	1	1	1	-0.19258164	-0.1467916	0.48180824	-0.10021618	1.2753218	⋯	-0.9481277	-0.8560700	1.1377171	-0.3378133	0.18457796	0.1073270	-0.15329655	-1.00827077	-0.1413894	MULTIPOLYGON (((526123.5 47...
2	17	26	10	1	1	-0.36188692	-0.3284478	1.32081073	0.44848018	-1.3865567	⋯	0.0000000	-0.2328068	-0.9881374	-0.4065854	-0.33477437	0.8522391	-0.26183757	0.04893220	-0.4444888	MULTIPOLYGON (((528744.7 46...
3	17	26	100	1	1	0.01083739	-1.4002199	0.23798535	-0.52201032	-1.3865567	⋯	0.0000000	-0.3822322	-1.5969460	-0.4065854	-0.22140412	1.0122262	-0.17316713	-0.05159845	-0.6491854	MULTIPOLYGON (((575382.6 46...
4	17	26	101	1	1	-1.36078806	-1.4002199	-0.22155305	0.07955613	-1.3865567	⋯	0.0000000	-1.5116073	-1.8689267	-0.4065854	0.02693072	1.3626742	0.02106336	-0.81142134	-0.9789365	MULTIPOLYGON (((540334.4 46...
5	17	26	102	1	1	0.10257030	-0.4581527	0.05129390	0.59299501	0.3880289	⋯	1.2637707	-0.3792937	0.4071804	-0.3116725	0.44354790	0.7414772	-0.09187444	0.80062025	0.4454417	MULTIPOLYGON (((521986.4 46...
6	17	26	102	1	2	0.64987009	0.6302157	-0.04065368	0.64469364	0.3880289	⋯	0.7211874	-0.3696338	0.7600790	-0.3203095	0.33117214	0.3141994	-0.30936550	0.44207919	0.2948932	MULTIPOLYGON (((519949.8 47...

Export

# CSV
write.csv(output_dataset, file.path(output_dir, "social_vulnerability_index_logrono_2021.csv"), row.names = FALSE)

# GeoJSON
st_write(output_dataset_geom, file.path(output_dir, "social_vulnerability_index_logrono_2021.geojson"), delete_dsn=TRUE)

# Shapefile
# Need to manually rename these fields, otherwise we get a shapefile creation error
names(output_dataset_geom)[names(output_dataset_geom) == 'early_childhood_boy'] <- 'erly_cld_b'
names(output_dataset_geom)[names(output_dataset_geom) == 'early_childhood_girl'] <- 'erly_cld_g'
names(output_dataset_geom)[names(output_dataset_geom) == 'age_middle_to_oldest_old_male'] <- 'age_old_m'
names(output_dataset_geom)[names(output_dataset_geom) == 'age_middle_to_oldest_old_female'] <- 'age_old_f'
st_write(output_dataset_geom, file.path(output_dir, "social_vulnerability_index_logrono_2021.shp"), append = FALSE)

Deleting source `../../3_outputs/Spain/Logrono/2021/social_vulnerability_index_logrono_2021.geojson' using driver `GeoJSON'
Writing layer `social_vulnerability_index_logrono_2021' to data source 
  `../../3_outputs/Spain/Logrono/2021/social_vulnerability_index_logrono_2021.geojson' using driver `GeoJSON'
Writing 343 features with 38 fields and geometry type Multi Polygon.

Warning message in abbreviate_shapefile_names(obj):
“Field names abbreviated for ESRI Shapefile driver”

Deleting layer `social_vulnerability_index_logrono_2021' using driver `ESRI Shapefile'
Writing layer `social_vulnerability_index_logrono_2021' to data source 
  `../../3_outputs/Spain/Logrono/2021/social_vulnerability_index_logrono_2021.shp' using driver `ESRI Shapefile'
Writing 343 features with 38 fields and geometry type Multi Polygon.

END