Composite accessibility index - Munster

# Load libraries
library(sf)
library(dplyr)
library(readr)

# Set working directory
setwd("C:/Users/eujoe/Documents/welbeing/Final_submission_2")
getwd()

[1] "C:/Users/eujoe/Documents/welbeing/Final_submission_2"

# List all files
list.files()

 [1] "betweenes_clean.csv"                   
 [2] "clossness_clean.csv"                   
 [3] "exercise.nb.html"                      
 [4] "exercise.Rmd"                          
 [5] "frequency_clean.csv"                   
 [6] "Grid.cpg"                              
 [7] "Grid.dbf"                              
 [8] "Grid.prj"                              
 [9] "Grid.qmd"                              
[10] "Grid.shp"                              
[11] "Grid.shx"                              
[12] "pedestrian_network_clean.csv"          
[13] "public_Transportation_access_clean.csv"

# Read the spatial grid
grid <- st_read("Grid.shp")

Reading layer `Grid' from data source 
  `C:\Users\eujoe\Documents\welbeing\Final_submission_2\Grid.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1331 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 395103.5 ymin: 5744177 xmax: 415704.7 ymax: 5768641
Projected CRS: ETRS89 / UTM zone 32N

betweenes <- read.csv(
  "betweenes_clean.csv",
  stringsAsFactors = FALSE,
  encoding = "UTF-8"
)

closeness <- read.csv(
  "clossness_clean.csv",
  stringsAsFactors = FALSE,
  fileEncoding = "UTF-8"
)

pedestrian <- read.csv(
  "pedestrian_network_clean.csv",
  stringsAsFactors = FALSE,
  fileEncoding = "UTF-8"
)

public_tr <- read.csv(
  "public_Transportation_access_clean.csv",
  stringsAsFactors = FALSE,
  fileEncoding = "UTF-8"
)

frequency <- read.csv(
  "frequency_clean.csv",
  stringsAsFactors = FALSE,
  fileEncoding = "UTF-8"
)

names(grid)

[1] "OBJECTID"   "GRID_ID"    "Shape_Leng" "Shape_Area" "geometry"  

names(betweenes)

[1] "GRID_ID"   "bet_value"

names(closeness)

[1] "GRID_ID"    "clos_value"

names(pedestrian)

[1] "GRID_ID"  "pd_value"

names(public_tr)

[1] "GRID_ID"   "pta_value"

names(frequency)

[1] "GRID_ID"   "frq_value"

# Check IDs are the same type
grid$GRID_ID        <- as.character(grid$GRID_ID)
betweenes$GRID_ID <- as.character(betweenes$GRID_ID)
closeness$GRID_ID   <- as.character(closeness$GRID_ID)
pedestrian$GRID_ID  <- as.character(pedestrian$GRID_ID)
public_tr$GRID_ID   <- as.character(public_tr$GRID_ID)
frequency$GRID_ID   <- as.character(frequency$GRID_ID)

# Join CSVs together
attributes <- betweenes %>%
  left_join(closeness,  by = "GRID_ID") %>%
  left_join(pedestrian, by = "GRID_ID") %>%
  left_join(public_tr,  by = "GRID_ID") %>%
  left_join(frequency,  by = "GRID_ID")

# Join to spatial grid
grid_joined <- grid %>%
  left_join(attributes, by = "GRID_ID")

# Drop geometry
grid_joined_df <- st_drop_geometry(grid_joined)

# Write to CSV
write.csv(
  grid_joined_df,
  "Grid_wellbeing.csv",
  row.names = FALSE
)

grid_joined_df <- st_drop_geometry(grid_joined)
summary(grid_joined$bet_value)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
       0        0       59  1099991   528882 45269913 

summary(grid_joined$clos_value)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
 0.00000  0.00000  0.06108  4.76393 10.25823 17.84093 

summary(grid_joined$frq_value)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    5.0   715.0   963.2   904.9  1113.9  1934.0     895 

sum(!is.finite(grid_joined$bet_value))

[1] 0

sum(!is.finite(grid_joined$clos_value))

[1] 0

sum(!is.finite(grid_joined$frq_value))

[1] 895

grid_joined <- grid_joined %>%
  mutate(
    bet_value  = ifelse(is.finite(bet_value),  bet_value,  0),
    clos_value = ifelse(is.finite(clos_value), clos_value, 0),
    frq_value = ifelse(is.finite(frq_value), frq_value, 0)
  )

library(spdep)

# Create neighbors
nb <- poly2nb(grid_joined, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Create lag variables INSIDE grid_joined
grid_joined <- grid_joined %>%
  mutate(
    bet_lag  = lag.listw(lw, bet_value,  zero.policy = TRUE),
    clos_lag = lag.listw(lw, clos_value, zero.policy = TRUE) ,
    frq_lag = lag.listw(lw, frq_value, zero.policy = TRUE)
    
    
  )

# Drop geometry (lags included)
df <- st_drop_geometry(grid_joined)

z_std <- function(x) {
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

df_std <- df %>%
  mutate(
    bet_z  = z_std(bet_lag),
    clos_z = z_std(clos_lag),
    pd_z   = z_std(pd_value),
    pta_z  = z_std(pta_value),
    frq_z  = z_std(frq_lag),
  )

names(df)

 [1] "OBJECTID"   "GRID_ID"    "Shape_Leng" "Shape_Area" "bet_value" 
 [6] "clos_value" "pd_value"   "pta_value"  "frq_value"  "bet_lag"   
[11] "clos_lag"   "frq_lag"   

summary(df$bet_lag)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 2.210e+01 1.077e+05 1.102e+06 1.625e+06 1.244e+07 

summary(df$clos_lag)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
 0.00000  0.02277  3.72398  4.78497  8.52862 16.13242 

summary(df$frq_lag)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.0   154.9   297.9   531.9  1308.7 

df_std <- df_std %>%
  mutate(
    accessibility_index =
      (bet_z + clos_z + pd_z + pta_z+frq_z) / 5
  )

write.csv(
  df_std,
  "Grid_accessibility_index_mixed.csv",
  row.names = FALSE
)

df_std <- df_std %>%
  mutate(
    accessibility_index_01 =
      (accessibility_index - min(accessibility_index, na.rm = TRUE)) /
      (max(accessibility_index, na.rm = TRUE) - min(accessibility_index, na.rm = TRUE))
  )

grid_access <- grid_joined %>%
  left_join(
    df_std %>% select(GRID_ID, accessibility_index_01),
    by = "GRID_ID"
  )

st_write(
  grid_access,
  "Grid_accessibility_index_mixed.gpkg",
  layer = "accessibility",
  delete_layer = TRUE
)

Deleting layer `accessibility' using driver `GPKG'
Writing layer `accessibility' to data source 
  `Grid_accessibility_index_mixed.gpkg' using driver `GPKG'
Writing 1331 features with 13 fields and geometry type Multi Polygon.

library(ggplot2)
ggplot(grid_access) +
  geom_sf(aes(fill = accessibility_index_01), color = NA) +
  scale_fill_viridis_c(
    name = "Accessibility\nIndex",
    option = "viridis"
  ) +
  labs(
    title = "Composite Accessibility Index",
    caption = "Higher values indicate better accessibility"
  ) +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank()
  )

library(tmap)
library(spdep)
library(dplyr)

tmap_mode("plot")

grid_access <- grid_joined %>%
  left_join(
    df_std %>%
      select(GRID_ID,
             bet_z, clos_z, pd_z, pta_z,
             accessibility_index,
             accessibility_index_01),
    by = "GRID_ID"
  )

# Z-score function
z_std <- function(x) {
  if (sd(x, na.rm = TRUE) == 0) {
    rep(0, length(x))
  } else {
    (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
  }
}

grid_access <- grid_access %>%
  mutate(
    bet_z      = z_std(bet_value),
    bet_z_lag  = z_std(bet_lag),
    clos_z     = z_std(clos_value),
    clos_z_lag = z_std(clos_lag),
    frq_z  = z_std(frq_value),
    frq_z_lag = z_std(frq_lag),
    
  )
names(grid_access)

 [1] "OBJECTID"               "GRID_ID"               
 [3] "Shape_Leng"             "Shape_Area"            
 [5] "bet_value"              "clos_value"            
 [7] "pd_value"               "pta_value"             
 [9] "frq_value"              "bet_lag"               
[11] "clos_lag"               "frq_lag"               
[13] "bet_z"                  "clos_z"                
[15] "pd_z"                   "pta_z"                 
[17] "accessibility_index"    "accessibility_index_01"
[19] "geometry"               "bet_z_lag"             
[21] "clos_z_lag"             "frq_z"                 
[23] "frq_z_lag"             

bet_z_map <- tm_shape(grid_access) +
  tm_polygons("bet_z", palette = "viridis", title = "Betweenness (z)") +
  tm_layout(frame = FALSE)

bet_z_lag_map <- tm_shape(grid_access) +
  tm_polygons("bet_z_lag", palette = "viridis", title = "Betweenness (lag + z)") +
  tm_layout(frame = FALSE)

clos_z_map <- tm_shape(grid_access) +
  tm_polygons("clos_z", palette = "viridis", title = "Closeness (z)") +
  tm_layout(frame = FALSE)

clos_z_lag_map <- tm_shape(grid_access) +
  tm_polygons("clos_z_lag", palette = "viridis", title = "Closeness (lag + z)") +
  tm_layout(frame = FALSE)

frq_z_map <- tm_shape(grid_access) +
  tm_polygons("frq_z", palette = "viridis", title = "Frequency (z)") +
  tm_layout(frame = FALSE)

frq_z_lag_map <- tm_shape(grid_access) +
  tm_polygons("frq_z_lag", palette = "viridis", title = "Frequency (lag + z)") +
  tm_layout(frame = FALSE)

tmap_arrange(
  bet_z_map, bet_z_lag_map,
  clos_z_map, clos_z_lag_map,
   frq_z_map, frq_z_lag_map,
  
  ncol = 2
)

tm_shape(grid_access) +
  tm_polygons(
    "pd_z",
    style = "quantile",
    n = 5,
    palette = "viridis",
    title = "Pedestrian Network (z)"
  ) +
  tm_layout(frame = FALSE)

tm_shape(grid_access) +
  tm_polygons(
    "pta_z",
    style = "quantile",
    n = 5,
    palette = "viridis",
    title = "Public Transport Access (z)"
  ) +
  tm_layout(frame = FALSE)

grid_access$accessibility_index_01[
  is.na(grid_access$accessibility_index_01)
] <- 0

nb <- poly2nb(grid_access, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

z <- scale(grid_access$accessibility_index_01)[, 1]
z[!is.finite(z)] <- 0

z_lag <- lag.listw(lw, z, zero.policy = TRUE)
z_lag[!is.finite(z_lag)] <- 0

lisa <- localmoran(
  grid_access$accessibility_index_01,
  lw,
  zero.policy = TRUE,
  na.action = na.exclude
)
grid_access$lisa_I <- lisa[, 1]
grid_access$lisa_p <- lisa[, ncol(lisa)]

# Local Moran's I
lisa <- localmoran(
  grid_access$accessibility_index_01,
  lw,
  zero.policy = TRUE,
  na.action = na.exclude
)

# Attach results safely
grid_access$lisa_I <- lisa[, 1]
grid_access$lisa_p <- lisa[, ncol(lisa)]

z <- scale(grid_access$accessibility_index_01)[,1]
z[!is.finite(z)] <- 0

z_lag <- lag.listw(lw, z, zero.policy = TRUE)
z_lag[!is.finite(z_lag)] <- 0

grid_access$lisa_cluster <- "Not significant"

grid_access$lisa_cluster[
  z > 0 & z_lag > 0 & grid_access$lisa_p <= 0.05
] <- "High–High"

grid_access$lisa_cluster[
  z < 0 & z_lag < 0 & grid_access$lisa_p <= 0.05
] <- "Low–Low"

grid_access$lisa_cluster[
  z > 0 & z_lag < 0 & grid_access$lisa_p <= 0.05
] <- "High–Low"

grid_access$lisa_cluster[
  z < 0 & z_lag > 0 & grid_access$lisa_p <= 0.05
] <- "Low–High"

library(tmap)

tm_shape(grid_access) +
  tm_polygons(
    "lisa_cluster",
    palette = c(
      "High–High" = "red",
      "Low–Low" = "blue",
      "High–Low" = "pink",
      "Low–High" = "lightblue",
      "Not significant" = "grey80"
    ),
    title = "LISA Clusters (p ≤ 0.05)"
  ) +
  tm_layout(
    frame = FALSE,
    legend.outside = TRUE
     )

tm_shape(grid_access) +
  tm_polygons(
    "accessibility_index_01",
  style = "fixed",
  breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1),
  palette = "viridis",
  title = "Accessibility Index"
)