diff --git a/.DS_Store b/.DS_Store
index 9dc14be..e322e24 100644
Binary files a/.DS_Store and b/.DS_Store differ
diff --git a/.Rproj.user/.DS_Store b/.Rproj.user/.DS_Store
index 843647e..6927363 100644
Binary files a/.Rproj.user/.DS_Store and b/.Rproj.user/.DS_Store differ
diff --git a/.Rproj.user/C57E093E/.DS_Store b/.Rproj.user/C57E093E/.DS_Store
index c0d7450..9f9803b 100644
Binary files a/.Rproj.user/C57E093E/.DS_Store and b/.Rproj.user/C57E093E/.DS_Store differ
diff --git a/.Rproj.user/C57E093E/pcs/files-pane.pper b/.Rproj.user/C57E093E/pcs/files-pane.pper
index e13fefa..0a62c77 100644
--- a/.Rproj.user/C57E093E/pcs/files-pane.pper
+++ b/.Rproj.user/C57E093E/pcs/files-pane.pper
@@ -5,5 +5,5 @@
             "ascending": true
         }
     ],
-    "path": "~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access"
+    "path": "~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/Greater_Manchester_Combined_Authority"
 }
\ No newline at end of file
diff --git a/.Rproj.user/C57E093E/pcs/source-pane.pper b/.Rproj.user/C57E093E/pcs/source-pane.pper
index 544e582..902cc6f 100644
--- a/.Rproj.user/C57E093E/pcs/source-pane.pper
+++ b/.Rproj.user/C57E093E/pcs/source-pane.pper
@@ -1,3 +1,3 @@
 {
-    "activeTab": 5
+    "activeTab": 0
 }
\ No newline at end of file
diff --git a/.Rproj.user/C57E093E/pcs/workbench-pane.pper b/.Rproj.user/C57E093E/pcs/workbench-pane.pper
index ab5e950..1893013 100644
--- a/.Rproj.user/C57E093E/pcs/workbench-pane.pper
+++ b/.Rproj.user/C57E093E/pcs/workbench-pane.pper
@@ -1,5 +1,5 @@
 {
     "TabSet1": 3,
-    "TabSet2": 0,
+    "TabSet2": 1,
     "TabZoom": {}
 }
\ No newline at end of file
diff --git a/.Rproj.user/C57E093E/sources/prop/6B25F063 b/.Rproj.user/C57E093E/sources/prop/6B25F063
index fe9deac..b18f217 100644
--- a/.Rproj.user/C57E093E/sources/prop/6B25F063
+++ b/.Rproj.user/C57E093E/sources/prop/6B25F063
@@ -1,6 +1,6 @@
 {
     "source_window_id": "",
     "Source": "Source",
-    "cursorPosition": "90,72",
-    "scrollLine": "77"
+    "cursorPosition": "96,0",
+    "scrollLine": "0"
 }
\ No newline at end of file
diff --git a/.Rproj.user/C57E093E/sources/prop/INDEX b/.Rproj.user/C57E093E/sources/prop/INDEX
index 8ffab5a..01b0fb9 100644
--- a/.Rproj.user/C57E093E/sources/prop/INDEX
+++ b/.Rproj.user/C57E093E/sources/prop/INDEX
@@ -26,7 +26,11 @@
 ~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FLiverpool_City_Region%2FInteractive_Map_LCR.R="869C046C"
 ~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F0.Decay_Params_WYCA.R="3DB85296"
 ~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F0.Rail_BODS_GTFS_files_WYCA.R="83AC253B"
+~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F1._Public_Transport_Job_Access_WYCA.R="66656B5E"
 ~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F1._TravelTimeMatrix_WYCA.R="626002FF"
+~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F2._Compile_Datasets_WYCA.R="F08A16E7"
+~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F3.a.Job_opportunities_WYCA.R="E6B7E610"
+~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F3.b.Population_Density_Map_WYCA.R="23B1CB72"
 ~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F3.c.PopulationDensityChange_WYCA.R="F5D9A97A"
 ~%2FLibrary%2FCloudStorage%2FGoogleDrive-sam.allwood3%40gmail.com%2FMy%20Drive%2FConsulting%2FUnemployment_Public_Transport_Access%2FWYCA%2F3.d.PTJA_WYCA.R="1A9FC226"
 ~%2FLibrary%2FCloudStorage%2FOneDrive-UniversityofGlasgow%2FConsulting%2FUnemployment_Public_Transport_Access%2F0.Decay_Params.R="F9521DCA"
diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths
index b83dabb..6a8665d 100644
--- a/.Rproj.user/shared/notebooks/paths
+++ b/.Rproj.user/shared/notebooks/paths
@@ -1,3 +1,11 @@
 /Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/.gitignore="C5705872"
 /Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/0.Decay_Params.R="13A251CE"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/0.Rail_BODS_GTFS_files.R="5B1F04D7"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/Greater_Manchester_Combined_Authority/1._Public_Transport_Job_Access.R="743BF62E"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/Liverpool_City_Region/1._TravelTimeMatrix_Liverpool.R="E3AEDDA5"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/WYCA/0.Decay_Params_WYCA.R="18908051"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/WYCA/1._Public_Transport_Job_Access_WYCA.R="FBC25968"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/WYCA/2._Compile_Datasets_WYCA.R="B0A0E905"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/WYCA/3.b.Population_Density_Map_WYCA.R="02CEBB0B"
+/Users/samallwood/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/WYCA/3.c.PopulationDensityChange_WYCA.R="172F1207"
 /Users/samallwood/Library/CloudStorage/OneDrive-UniversityofGlasgow/Consulting/Unemployment_Public_Transport_Access/Greater_Manchester_Combined_Authority/3.b.Population_Density_Map.R="24D36CE0"
diff --git a/0.Rail_BODS_GTFS_files.R b/0.Rail_BODS_GTFS_files.R
index 28d142b..ef9ff4b 100644
--- a/0.Rail_BODS_GTFS_files.R
+++ b/0.Rail_BODS_GTFS_files.R
@@ -12,8 +12,8 @@ library(here)
 # 2. Load and Filter Datasets ------------------------------------------------------------
 # Combined Authority Boundary + buffer
 Boundaries <- read_sf("../Data/CAUTH_DEC_2023_EN_BFC.shp")
-CA_boundary <- Boundaries %>% filter(CAUTH23NM == "Greater Manchester" | 
-                                       CAUTH23NM =="Liverpool City Region" | 
+CA_boundary <- Boundaries %>% filter(#CAUTH23NM == "Greater Manchester" | 
+                                     #  CAUTH23NM =="Liverpool City Region" | 
                                        CAUTH23NM =="West Yorkshire") %>%
   st_transform(4326) 
 #mapview(CA_boundary)
@@ -30,7 +30,7 @@ combined_area <- st_union(buffered_CA_boundary)
 BODS <- gtfstools::read_gtfs("../Data/r5r_data/itm_all_gtfs.zip") # Download from BODS
 # Filter BODS stops in CAs
 BODS_CA <- filter_by_spatial_extent(BODS, buffered_CA_boundary)
-gtfstools::write_gtfs(BODS_CA, "../Data/r5r_data/BODS_CA.gtfs.zip")
+gtfstools::write_gtfs(BODS_CA, "../Data/r5r_data/CA/BODS_CA.gtfs.zip")
 
 # Transform ATOC (rail data) ----------------------------------------------
 # Uses UK2GTFS V. ‘0.1.1’
@@ -61,10 +61,11 @@ na_columns <- colSums(is.na(stops)) > 0
 print(na_columns)
 
 # Update the GTFS stops data
-ttis272$stops <- stops
+ttis272_corrected <- ttis272
+ttis272_corrected$stops <- stops
 
 ## Force valid. This function does not fix problems, it just removes them
-ttis272_gtfs <- UK2GTFS::gtfs_force_valid(ttis272)
+ttis272_gtfs <- UK2GTFS::gtfs_force_valid(ttis272_corrected)
 UK2GTFS::gtfs_validate_internal(ttis272_gtfs)
 
 ## Compare original and valid
@@ -80,6 +81,8 @@ count(gtfs_diff, stop_id)
 # Trips affected
 unique(gtfs_diff$trip_id)
 
+# map the eliminated stop_ids somehow
+
 UK2GTFS::gtfs_write(ttis272_gtfs, 
                     quote = TRUE, 
                     folder = "../Data/r5r_data", 
@@ -88,7 +91,7 @@ rail_gtfstools <- gtfstools::read_gtfs("../Data/r5r_data/rail.gtfs.zip")
 rail_CA_gtfs <- filter_by_spatial_extent(rail_gtfstools, buffered_CA_boundary)
 UK2GTFS::gtfs_validate_internal(rail_CA_gtfs)
 
-gtfstools::write_gtfs(rail_CA_gtfs, "../Data/r5r_data/rail_CA.gtfs.zip")
+gtfstools::write_gtfs(rail_CA_gtfs, "../Data/r5r_data/CA/rail_CA.gtfs.zip")
 
 
 latest_validator <- download_validator(tempdir())
diff --git a/Liverpool_City_Region/.DS_Store b/Liverpool_City_Region/.DS_Store
index fecd53f..0b183b1 100644
Binary files a/Liverpool_City_Region/.DS_Store and b/Liverpool_City_Region/.DS_Store differ
diff --git a/WYCA/.DS_Store b/WYCA/.DS_Store
index 94b0e95..b2e2be8 100644
Binary files a/WYCA/.DS_Store and b/WYCA/.DS_Store differ
diff --git a/WYCA/1._Public_Transport_Job_Access_WYCA.R b/WYCA/1._Public_Transport_Job_Access_WYCA.R
index 45fcc44..94f0cc7 100644
--- a/WYCA/1._Public_Transport_Job_Access_WYCA.R
+++ b/WYCA/1._Public_Transport_Job_Access_WYCA.R
@@ -16,10 +16,8 @@ options(timeout = 1000)
 options(scipen=999)
 
 # 2. Load and Filter Datasets ------------------------------------------------------------
-setwd("~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/WYCA")
-
 # WYCA Boundary + buffer
-Boundaries <- read_sf("../../Data/CAUTH_DEC_2023_EN_BFC.shp")
+Boundaries <- read_sf("../Data/CAUTH_DEC_2023_EN_BFC.shp")
 WYCA_boundary <- Boundaries %>% filter(CAUTH23NM == "West Yorkshire") %>%
   st_transform(4326) 
 WYCA_bound_small_buffer <- WYCA_boundary %>% st_buffer(dist=25)
@@ -29,7 +27,7 @@ buffered_WYCA_boundary <- st_buffer(WYCA_boundary, dist = 20000) %>%
   st_transform(4326)
 
 # Read Local Authority District (LAD) boundaries
-LADs <- read_sf("../../Data/LAD_DEC_2021_GB_BFC.shp") %>%
+LADs <- read_sf("../Data/LAD_DEC_2021_GB_BFC.shp") %>%
   st_transform(4326)
 # Filter LADs within WYCA
 LADs_WYCA <- LADs %>% filter(as.vector(st_within(., WYCA_bound_small_buffer, sparse = FALSE))) %>% 
@@ -37,7 +35,7 @@ LADs_WYCA <- LADs %>% filter(as.vector(st_within(., WYCA_bound_small_buffer, spa
 # mapview(LADs_WYCA)+mapview(WYCA_boundary)
 
 # Read LSOA pop-weighted centroids
-lsoas <- st_read("../../Data/LSOA_PopCentroids_EW_2021_V3.shp") %>%
+lsoas <- st_read("../Data/LSOA_PopCentroids_EW_2021_V3.shp") %>%
   st_transform(4326)
 # Filter LSOA PW-centroids
 lsoas <- st_make_valid(lsoas)
@@ -52,7 +50,7 @@ lsoa_PWC_within_WYCA_buffer$geometry <- st_transform(lsoa_PWC_within_WYCA_buffer
 lsoa_PWC_within_WYCA$geometry <- st_transform(lsoa_PWC_within_WYCA$geometry, 4326)
 
 #Read LSOA boundaries
-lsoa_boundaries <- st_read("../../Data/LSOA_2021_EW_BFC_V8.shp") %>%
+lsoa_boundaries <- st_read("../Data/LSOA_2021_EW_BFC_V8.shp") %>%
   st_transform(4326)
 lsoa_boundaries$geometry <- st_make_valid(lsoa_boundaries$geometry)
 
@@ -69,17 +67,17 @@ lsoa_boundaries_within_WYCA_buffer$LSOA_area_km2 <- lsoa_boundaries_within_WYCA_
 lsoa_boundaries_within_WYCA$LSOA_area <- st_area(lsoa_boundaries_within_WYCA$geometry)
 lsoa_boundaries_within_WYCA$LSOA_area_km2 <- lsoa_boundaries_within_WYCA$LSOA_area/1000000
 # Read Towns and City boundaries
-towns <- st_read("../../Data/TCITY_2015_EW_BGG_V2.shp") %>%
+towns <- st_read("../Data/TCITY_2015_EW_BGG_V2.shp") %>%
   st_transform(4326) 
 towns$geometry <- st_make_valid(towns$geometry)
 towns_within_WYCA_buffer <- towns[st_within(towns, buffered_WYCA_boundary, sparse = FALSE), ]
 towns_within_WYCA <- towns[st_within(towns, WYCA_bound_small_buffer, sparse = FALSE), ]
 
 towns_manual_WYCA <- data.frame( 
-  id = c("Barnsley", "Bradford", "Burnley", "Bury", "Doncaster", "Halifax", "Harrogate", "Huddersfield", "Leeds", "Oldham"),
-  lon = c(-1.479237, -1.756356, -2.245021, -2.297645, -1.131974, -1.859140, -1.537945, -1.782488, -1.545146, -2.112658), 
-  lat = c(53.553464, 53.794762, 53.789720, 53.592658, 53.524176, 53.723237,53.996504, 53.646861, 53.796904, 53.540754),
-  buffer = c("Y", "N", "Y", "Y", "Y", "N", "Y", "N", "N", "Y")
+  id = c("Barnsley", "Bradford", "Burnley", "Bury", "Doncaster", "Halifax", "Harrogate", "Huddersfield", "Leeds", "Oldham", "Rochdale", "Rotherham", "Wakefield", "York"),
+  lon = c(-1.479237, -1.756356, -2.245021, -2.297645, -1.131974, -1.859140, -1.537945, -1.782488, -1.545146, -2.112658, -2.148727, -1.355379, -1.505890, -1.080324),
+  lat = c(53.553464, 53.794762, 53.789720, 53.592658, 53.524176, 53.723237,53.996504, 53.646861, 53.796904, 53.540754, 53.612109, 53.432603, 53.676009, 53.959055),
+  buffer = c("Y", "N", "Y", "Y", "Y", "N", "Y", "N", "N", "Y", "Y", "Y", "N", "Y")
 ) 
 
 # Convert the dataframe to an sf object, assuming WGS 84 (EPSG: 4326) coordinate reference system
@@ -117,7 +115,7 @@ lsoa_with_town_info_min <- lsoa_with_town_info %>%  select(LSOA21CD,
 #mapview(towns_within_WYCA_buffer)+mapview(WYCA_boundary)+mapview(towns_centroids)+mapview(towns_manual_WYCA)
 #mapview(lsoa_with_town_info)+ mapview(towns_with_buffer)
 
-# Liverpool City Centre - LOCATION not super-important as not used for any calculation. Arbitrary city centre point
+# Leeds City Centre - LOCATION not super-important as not used for any calculation. Arbitrary city centre point
 Leeds_CC <- 
   data.frame(id = "Leeds_CC", lat = 53.796904, lon = -1.545146) %>% 
   st_as_sf(coords = c("lon", "lat"), crs = 4326)
@@ -125,7 +123,7 @@ Leeds_CC_point <-
   geom_sf(data = Leeds_CC, shape = 21, fill = 'white', size = 2)
 mapview(Leeds_CC)
 # Employed Residents
-Employment <-read_csv("../../Data/census2021-ts066-lsoa.csv") %>%
+Employment <-read_csv("../Data/census2021-ts066-lsoa.csv") %>%
   dplyr::select("geography",
                 "geography code",
                 "Economic activity status: Economically active (excluding full-time students)") %>%
@@ -134,7 +132,7 @@ Employment <-read_csv("../../Data/census2021-ts066-lsoa.csv") %>%
          "Economically_active" = "Economic activity status: Economically active (excluding full-time students)")
 
 
-# Load and validate BODS GTFS data - for diagnostics only
+# Load and validate BODS GTFS data - code for diagnostics only
 #BODS_WYCA <- read_gtfs("Data/GTFS_Data/r5r/BODS_WYCA.gtfs.zip")
 #BODS_WYCA_valid <- tidytransit::validate_gtfs(BODS_WYCA, warnings = TRUE)
 #view(BODS_manch_valid)
@@ -152,7 +150,7 @@ Employment <-read_csv("../../Data/census2021-ts066-lsoa.csv") %>%
 # 4. R5R Network Setup -------------------------------------------------------
 # Get input files - will load BODS_MANCH and BODS_WYCA - consider simplifying if
 # too slow.
-input_files <- list.files("../../Data/r5r_data", 
+input_files <- list.files("../Data/r5r_data/CA", 
   recursive = TRUE,
   pattern = 'gtfs\\.zip$|pbf$',
   full.names = TRUE)
@@ -163,7 +161,7 @@ input_files
 # Read multi-modal network
 gc() #this setup step requires quite a bit of memory, so best to gc first
 options(java.parameters = "-Xmx8G")
-r5r_WYCA <- setup_r5(data_path = "../../Data/r5r_data", 
+r5r_WYCA <- setup_r5(data_path = "../Data/r5r_data/CA", 
                      verbose=TRUE,
                      overwrite = FALSE) # make sure you leave this to run after it indicates it's finished
 
@@ -182,17 +180,18 @@ print(r5r_WYCA)
 departure_datetime <- as.POSIXct("2024-05-21 08:00:00") 
 mode = c( "TRANSIT") # note function always includes 'WALK' anyway
 walk_speed = 4.32
-max_duration = 2000L
+max_duration = 20000L
 max_walk_time = 30L
 
 WYCA_TT_CC <- 
   travel_time_matrix(
     r5r_core = r5r_WYCA, 
     origins = lsoa_PWC_within_WYCA, 
-    destinations = WYCA_CC, 
+    destinations = Leeds_CC, 
     mode = mode,
     departure_datetime = departure_datetime,
     walk_speed = walk_speed,
+    max_trip_duration = max_duration,
     max_walk_time = max_walk_time,
     verbose = FALSE,
     progress = TRUE) %>%
@@ -204,8 +203,9 @@ WYCA_Emp_ettm <-
   expanded_travel_time_matrix(
     r5r_core = r5r_WYCA, 
     origins = lsoa_PWC_within_WYCA, 
-    destinations = WYCA_CC, 
+    destinations = Leeds_CC, 
     mode = mode,
+    max_trip_duration = max_duration,
     departure_datetime = departure_datetime,
     walk_speed = walk_speed,
     max_trip_duration = 50,
@@ -222,7 +222,7 @@ WYCA_Emp_ttm <-
     mode = mode,
     departure_datetime = departure_datetime,
     walk_speed = walk_speed,
-    max_trip_duration = 500,
+    max_trip_duration = max_duration,
     max_walk_time = max_walk_time,
     verbose = FALSE,
     progress = TRUE) 
@@ -286,13 +286,13 @@ WYCA_ttm_walk <-
 
 # BRES --------------------------------------------------------------------
 # Load employment data from BRES. Note BRES uses 2011 LSOA codes, so lookup for 2021 conversion required
-BRES <- read_csv("../../Data/Business_reg_Emp_Surv(BRES)2021.csv", skip = 8) %>% 
+BRES <- read_csv("../Data/Business_reg_Emp_Surv(BRES)2021.csv", skip = 8) %>% 
   separate("...1", into = c("LSOA11CD", "LSOA11NM"), sep = " : ") %>%
   rename("Employed" = "...2") %>%
   dplyr::select(-c("...3", "...4","...5")) %>%
   na.omit()
 # Load 2011-2021 LSOA lookup table
-LSOA_lookup <- read_csv("../../Data/LSOA_(2011)_to_LSOA_(2021)_to_Local_Authority_District_(2022)_Lookup_for_England_and_Wales.csv")
+LSOA_lookup <- read_csv("../Data/LSOA_(2011)_to_LSOA_(2021)_to_Local_Authority_District_(2022)_Lookup_for_England_and_Wales.csv")
 
 # Lookup correct LSOA codes
 # left join will duplicate 2011 employment numbers if LSOA is merged in 2021. 
@@ -314,7 +314,7 @@ job_access <- gravity(
   travel_cost="travel_time_p50",
   land_use_data=BRES_2021_corrected,
   opportunity="Employed",
-  decay_function=decay_logistic(39.7,12.6), #note mean=39.7 and sd=12.6 are calculated in Decay_Params.R file.
+  decay_function=decay_logistic(40.3, 12.8), #note mean=40.3 and sd=12.8 are calculated in 0.Decay_Params_WYCA.R file.
   fill_missing_ids = TRUE) %>%
   rename("PT_Job_Access_Index" = "Employed")
 
@@ -323,7 +323,7 @@ job_access_bus <- gravity(
   travel_cost="travel_time_p50",
   land_use_data=BRES_2021_corrected,
   opportunity="Employed",
-  decay_function=decay_logistic(39.7,12.6), #note mean=39.7 and sd=12.6 are calculated in Decay_Params.R file.
+  decay_function=decay_logistic(40.3, 12.8), #note mean=40.3 and sd=12.8 are calculated in 0.Decay_Params_WYCA.R file.
   fill_missing_ids = TRUE) %>%
   rename("PT_Job_Access_Index_Bus" = "Employed")
 
@@ -332,7 +332,7 @@ job_access_train <- gravity(
   travel_cost="travel_time_p50",
   land_use_data=BRES_2021_corrected,
   opportunity="Employed",
-  decay_function=decay_logistic(39.7,12.6), #note mean=39.7 and sd=12.6 are calculated in Decay_Params.R file.
+  decay_function=decay_logistic(40.3, 12.8), #note mean=40.3 and sd=12.8 are calculated in 0.Decay_Params_WYCA.R file.
   fill_missing_ids = TRUE) %>%
   rename("PT_Job_Access_Index_Train" = "Employed")
 
@@ -341,7 +341,7 @@ job_access_walk <- gravity(
   travel_cost="travel_time_p50",
   land_use_data=BRES_2021_corrected,
   opportunity="Employed",
-  decay_function=decay_logistic(39.7,12.6), #note mean=39.7 and sd=12.6 are calculated in Decay_Params.R file.
+  decay_function=decay_logistic(40.3, 12.8), #note mean=40.3 and sd=12.8 are calculated in 0.Decay_Params_WYCA.R file.
   fill_missing_ids = TRUE) %>%
   rename("PT_Job_Access_Index_Walk" = "Employed")
 
@@ -359,7 +359,7 @@ opp <- data.frame(
   travel_cost="travel_time_p50",
   land_use_data=opp,
   opportunity="Employed",
-  decay_function=decay_logistic(39.7,12.6), #note mean=39.7 and sd=12.6 are calculated in Decay_Params.R file.
+  decay_function=decay_logistic(40.3, 12.8), #note mean=40.3 and sd=12.8 are calculated in 0.Decay_Params_WYCA.R file.
   fill_missing_ids = TRUE) %>%
   rename("PT_Job_Access_Index" = "Employed"))
 
@@ -384,7 +384,7 @@ demand_potential <- gravity(
   travel_cost="travel_time_p50",
   land_use_data=Employment,
   opportunity="Economically_active",
-  decay_function=decay_logistic(39.7,12.6), #note mean=39.7 and sd=12.6 are calculated in Decay_Params.R file.
+  decay_function=decay_logistic(40.3, 12.8), #note mean=40.3 and sd=12.8 are calculated in 0.Decay_Params_WYCA.R file.
   fill_missing_ids = TRUE) %>%
   rename("Demand_potential" = "Economically_active") %>%
   left_join(BRES_2021_corrected, by = "id") %>%
@@ -397,7 +397,7 @@ job_access_demand <- gravity(
   travel_cost="travel_time_p50",
   land_use_data=demand_potential,
   opportunity="Jobs_over_demand_potential",
-  decay_function=decay_logistic(39.7,12.6), #note mean=39.7 and sd=12.6 are calculated in Decay_Params.R file.
+  decay_function=decay_logistic(40.3, 12.8), #note mean=40.3 and sd=12.8 are calculated in 0.Decay_Params_WYCA.R file.
   fill_missing_ids = TRUE) %>%
   rename("PT_Job_Access_Index_demand" = "Jobs_over_demand_potential")
 
@@ -414,5 +414,5 @@ WYCA_dataset <- WYCA_TT_CC %>%
   left_join(job_access_demand, by = c("LSOA21CD"="id")) 
 
 # Write dataset to shapefile
-st_write(WYCA_dataset, "../../Data/WYCA_dataset.shp", append=FALSE)
-st_write(towns_manual_WYCA, "../../Data/towns_centres_WYCA.shp", append=FALSE)
+st_write(WYCA_dataset, "../Data/WYCA_dataset.shp", append=FALSE)
+st_write(towns_manual_WYCA, "../Data/towns_centres_WYCA.shp", append=FALSE)
diff --git a/WYCA/2._Compile_Datasets_WYCA.R b/WYCA/2._Compile_Datasets_WYCA.R
index 485d39c..dc58239 100644
--- a/WYCA/2._Compile_Datasets_WYCA.R
+++ b/WYCA/2._Compile_Datasets_WYCA.R
@@ -11,12 +11,9 @@ library(classInt)
 library(readODS)
 library(accessibility)
 
-setwd("~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/Liverpool_City_Region")
-
 # 2. Load Data ------------------------------------------------------------
-
-# Load Liverpool Geo data from file
-LCR_dataset <- read_sf("../../Data/LCR_dataset.shp") %>%
+# Load WYCA Geo data from file
+WYCA_dataset <- read_sf("../Data/WYCA_dataset.shp") %>%
                                    rename("LSOA_Code" = "LSOA21C",
                                           "LSOA_Name" = "LSOA21N",
                                           "PT_Job_Access_Index" = "PT_Jb_A_I",
@@ -28,41 +25,41 @@ LCR_dataset <- read_sf("../../Data/LCR_dataset.shp") %>%
                                           "Traveltime_empcent" = "Tr__E_C",
                                           "Traveltime_CC" = "Trvl_CC",
                                           "PT_Job_Access_Index_Demand" = "PT_Jb_A_I_")
-towns_centres_LCR <- read_sf("../../Data/towns_centres_LCR.shp") 
+towns_centres_WYCA <- read_sf("../Data/towns_centres_WYCA.shp") 
 
-# LCR Boundary + buffer
-Boundaries <- read_sf("../../Data/CAUTH_DEC_2023_EN_BFC.shp")
-LCR_boundary <- Boundaries %>% filter(CAUTH23NM == "Liverpool City Region") %>%
+# WYCA Boundary + buffer
+Boundaries <- read_sf("../Data/CAUTH_DEC_2023_EN_BFC.shp")
+WYCA_boundary <- Boundaries %>% filter(CAUTH23NM == "West Yorkshire") %>%
   st_transform(4326) 
-LCR_bound_small_buffer <- LCR_boundary %>% st_buffer(dist=25)
+WYCA_bound_small_buffer <- WYCA_boundary %>% st_buffer(dist=25)
 
-towns_centres_LCR <- towns_centres_LCR %>% filter(as.vector(st_within(., LCR_bound_small_buffer, sparse = FALSE))) %>% 
+towns_centres_WYCA <- towns_centres_WYCA %>% filter(as.vector(st_within(., WYCA_bound_small_buffer, sparse = FALSE))) %>% 
   st_transform(4326)
 # Read Local Authority District (LAD) boundaries
-LADs <- read_sf("../../Data/LAD_DEC_2021_GB_BFC.shp") %>%
+LADs <- read_sf("../Data/LAD_DEC_2021_GB_BFC.shp") %>%
   st_transform(4326)
-# Filter LADs within LCR
-LADs_LCR <- LADs %>% filter(as.vector(st_within(., LCR_bound_small_buffer, sparse = FALSE))) %>% 
+# Filter LADs within WYCA
+LADs_WYCA <- LADs %>% filter(as.vector(st_within(., WYCA_bound_small_buffer, sparse = FALSE))) %>% 
   st_transform(4326)
 
 # Output area lookup
-Output_Area_Lookup <- read_csv("../../Data/Output_Area_to_Lower_layer_Super_Output_Area_to_Middle_layer_Super_Output_Area_to_Local_Authority_District_(December_2021)_Lookup_in_England_and_Wales_v3.csv") %>%
+Output_Area_Lookup <- read_csv("../Data/Output_Area_to_Lower_layer_Super_Output_Area_to_Middle_layer_Super_Output_Area_to_Local_Authority_District_(December_2021)_Lookup_in_England_and_Wales_v3.csv") %>%
   dplyr::select(LSOA21CD, MSOA21CD, LAD22CD) %>%
   distinct(LSOA21CD, .keep_all=TRUE)
 
 # Jobcentres Data - downloaded and filtered in TravelTimeMatrix.R
-Jobcentre_plus_geocoded <- read_csv("../../Data/Jobcentre_locations_geocoded.csv")
+Jobcentre_plus_geocoded <- read_csv("../Data/Jobcentre_locations_geocoded.csv")
 # Convert job centres to simple feature 
 Jobcentre_plus_geocoded <- Jobcentre_plus_geocoded %>% 
   st_as_sf(coords = c("lon", "lat"), crs = 4326)
-# Filter jobcentres within LCR
-Jobcentres_within_LCR <- Jobcentre_plus_geocoded %>% 
+# Filter jobcentres within WYCA
+Jobcentres_within_WYCA <- Jobcentre_plus_geocoded %>% 
   st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
-  filter(as.vector(st_within(., LCR_boundary, sparse = FALSE)))
+  filter(as.vector(st_within(., WYCA_boundary, sparse = FALSE)))
 
 
 # Unemployment Levels
-Employment <-read_csv("../../Data/census2021-ts066-lsoa.csv") %>%
+Employment <-read_csv("../Data/census2021-ts066-lsoa.csv") %>%
   dplyr::select("geography",
          "geography code",
          "Economic activity status: Total: All usual residents aged 16 years and over",
@@ -77,7 +74,7 @@ Employment <-read_csv("../../Data/census2021-ts066-lsoa.csv") %>%
 
 
 #Car ownership rate
-Car_ownership <- read_csv("../../Data/census2021-ts045-lsoa.csv") %>%
+Car_ownership <- read_csv("../Data/census2021-ts045-lsoa.csv") %>%
   dplyr::select("geography",
          "geography code",
          "Number of cars or vans: Total: All households",
@@ -89,7 +86,7 @@ Car_ownership <- read_csv("../../Data/census2021-ts045-lsoa.csv") %>%
   mutate(No_car_rate = No_Cars*100/Total_Households)
 
 # Ethnicity
-Ethnicity <- read_csv("../../Data/census2021-ts021-lsoa.csv", skip=6) %>%
+Ethnicity <- read_csv("../Data/census2021-ts021-lsoa.csv", skip=6) %>%
   dplyr::select("2021 super output area - lower layer",
          "White",
          "%...11") %>%
@@ -100,14 +97,14 @@ Ethnicity <- read_csv("../../Data/census2021-ts021-lsoa.csv", skip=6) %>%
   mutate(LSOA_Code = sapply(strsplit(as.character(LSOA_Code), " : "), "[", 1))
 
 # Single-parent household
-SingleParent <- read_csv("../../Data/census2021-ts003-lsoa.csv", skip=6) %>%
+SingleParent <- read_csv("../Data/census2021-ts003-lsoa.csv", skip=6) %>%
   rename("LSOA_Code" = "2021 super output area - lower layer",
          "Single_Parent_Households" = "2021") %>%
   drop_na() %>%
   mutate(LSOA_Code = sapply(strsplit(as.character(LSOA_Code), " : "), "[", 1))
 
 # NS-Socio-Economic Classification
-NS_SEC <- read_csv("../../Data/census2021-ts062-lsoa.csv", skip=0) %>%
+NS_SEC <- read_csv("../Data/census2021-ts062-lsoa.csv", skip=0) %>%
   rename(LSOA_Code = "geography code",
          Total = "National Statistics Socio-economic Classification (NS-SEC): Total: All usual residents aged 16 years and over",
          L1_L2_L3 = "National Statistics Socio-economic Classification (NS-SEC): L1, L2 and L3 Higher managerial, administrative and professional occupations",
@@ -136,7 +133,7 @@ NS_SEC <- read_csv("../../Data/census2021-ts062-lsoa.csv", skip=0) %>%
 
 
 # Qualification levels
-Qualifications <- read_csv("../../Data/census2021-ts067-lsoa.csv") %>%
+Qualifications <- read_csv("../Data/census2021-ts067-lsoa.csv") %>%
   dplyr::select("geography code",
          "Highest level of qualification: Total: All usual residents aged 16 years and over",
          "Highest level of qualification: No qualifications",
@@ -158,7 +155,7 @@ Qualifications <- read_csv("../../Data/census2021-ts067-lsoa.csv") %>%
            Low_Qualified/Total_residents_over_16*100)
 
 # Urban rural classifications
-Urban_Rural <- read_ods("../../Data/Rural_Urban_Classification_2011_lookup_tables_for_small_area_geographies.ods", 
+Urban_Rural <- read_ods("../Data/Rural_Urban_Classification_2011_lookup_tables_for_small_area_geographies.ods", 
                         sheet = "LSOA11", 
                         skip=2) %>%
   rename("LSOA11CD" = "Lower Super Output Area 2011 Code",
@@ -169,7 +166,7 @@ Urban_Rural <- read_ods("../../Data/Rural_Urban_Classification_2011_lookup_table
           "Urban_Rural") 
 
 # Load 2011-2021 LSOA lookup table
-LSOA_lookup <- read_csv("../../Data/LSOA_(2011)_to_LSOA_(2021)_to_Local_Authority_District_(2022)_Lookup_for_England_and_Wales.csv")
+LSOA_lookup <- read_csv("../Data/LSOA_(2011)_to_LSOA_(2021)_to_Local_Authority_District_(2022)_Lookup_for_England_and_Wales.csv")
 # Lookup 2021 LSOA codes
 # left join 2011 Urban Rural classification if LSOA is split into two in 2021. 
 # Split LSOAs - population split between new LSOAs generated in 2011
@@ -182,7 +179,7 @@ Urban_Rural_2021 <- Urban_Rural_2021 %>%
 
 # 3. Data Cleaning --------------------------------------------------------
 # Join datasets
-LCR_dataset_full <- left_join(LCR_dataset, Employment, by = c("LSOA_Code" = "id"), keep=FALSE) %>%
+WYCA_dataset_full <- left_join(WYCA_dataset, Employment, by = c("LSOA_Code" = "id"), keep=FALSE) %>%
   left_join(Car_ownership, by = c("LSOA_Code" = "id"), keep=FALSE) %>%
   left_join(Ethnicity, by = c("LSOA_Code" = "LSOA_Code"), keep=FALSE) %>%
   left_join(SingleParent, by = c("LSOA_Code" = "LSOA_Code"), keep=FALSE) %>%
@@ -196,9 +193,9 @@ LCR_dataset_full <- left_join(LCR_dataset, Employment, by = c("LSOA_Code" = "id"
   rename(LAD22CD = LAD22CD.x)
 
 # Write dataset to csv
-write_csv(LCR_dataset_full, "../../Data/LCR_dataset_full.csv")
+write_csv(WYCA_dataset_full, "../Data/WYCA_dataset_full.csv")
 
-LCR_dataset_full %>% dplyr::select ( -c(
+WYCA_dataset_full %>% dplyr::select ( -c(
                                    "Population (16 and over)",
                                    Total_Households,
                                    No_Cars, 
@@ -213,18 +210,18 @@ LCR_dataset_full %>% dplyr::select ( -c(
                                    Level_2_qual,
                                    Apprent_qual
                                 )) %>%
-  st_write( "../../Data/LCR_dataset_full_sf.shp", append = FALSE)
+  st_write( "../Data/WYCA_dataset_full_sf.shp", append = FALSE)
 
 # Calculate higher-level unemployment rates
 Overall_unemp_rate <- sum(Employment$Unemployed) / sum(Employment$"Economically_active") *100
-LCR_unemp_rate <- sum(LCR_dataset_full$Unemployed) / sum(LCR_dataset_full$"Economically_active") *100
-Town_unemp_rate <- LCR_dataset_full %>%
+WYCA_unemp_rate <- sum(WYCA_dataset_full$Unemployed) / sum(WYCA_dataset_full$"Economically_active") *100
+Town_unemp_rate <- WYCA_dataset_full %>%
   filter(twnsbrb == "Urban") %>%
   summarise(Unemployed_sum = sum(Unemployed),
             Economically_active_sum = sum(Economically_active)) %>%
   mutate(Unemployment_rate = (Unemployed_sum / Economically_active_sum) * 100) %>%
   pull(Unemployment_rate)
-Suburb_unemp_rate <- LCR_dataset_full %>%
+Suburb_unemp_rate <- WYCA_dataset_full %>%
   filter(twnsbrb == "Suburb") %>%
   summarise(Unemployed_sum = sum(Unemployed),
             Economically_active_sum = sum(Economically_active)) %>%
diff --git a/WYCA/3.a.Job_opportunities_WYCA.R b/WYCA/3.a.Job_opportunities_WYCA.R
index e2ed164..f4358cf 100644
--- a/WYCA/3.a.Job_opportunities_WYCA.R
+++ b/WYCA/3.a.Job_opportunities_WYCA.R
@@ -1,3 +1,4 @@
+## Script to create a map of job opportunities in WYCA
 # 1. Setup ----------------------------------------------------------------
 library(tidyverse)
 library(knitr)
@@ -9,12 +10,9 @@ library(ggplot2)
 library(classInt)
 library(ggspatial)
 
-setwd("~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/Liverpool_City_Region")
-
 # 2. Load Data ------------------------------------------------------------
-
-# Load Liverpool Geo data from file
-LCR_dataset_jobs <- read_sf("../../Data/LCR_dataset.shp") %>%
+# Load WYCA Geo data from file
+WYCA_dataset_jobs <- read_sf("../Data/WYCA_dataset.shp") %>%
   rename("LSOA_Code" = "LSOA21C",
          "LSOA_Name" = "LSOA21N",
          "PT_Job_Access_Index" = "PT_Jb_A_I",
@@ -22,48 +20,56 @@ LCR_dataset_jobs <- read_sf("../../Data/LCR_dataset.shp") %>%
          "LSOA_Area" = "LSOA__2",
          "Traveltime_empcent" = "Tr__E_C")
 
-# LCR Boundary + buffer
-Boundaries <- read_sf("../../Data/CAUTH_DEC_2023_EN_BFC.shp")
-LCR_boundary <- Boundaries %>% filter(CAUTH23NM == "Liverpool City Region") %>%
+# WYCA Boundary + buffer
+Boundaries <- read_sf("../Data/CAUTH_DEC_2023_EN_BFC.shp")
+WYCA_boundary <- Boundaries %>% filter(CAUTH23NM == "West Yorkshire") %>%
   st_transform(4326) 
-LCR_bound_small_buffer <- LCR_boundary %>% st_buffer(dist=1000)
+WYCA_bound_small_buffer <- WYCA_boundary %>% st_buffer(dist=1000)
 
 # Town centres
-towns_centres_LCR <- read_sf("../../Data/towns_centres_LCR.shp") %>% 
+towns_centres_WYCA <- read_sf("../Data/towns_centres_WYCA.shp") %>% 
   st_transform(4326) %>%
   filter(buffer == "N") 
 
-#Town boundaries
-towns <- st_read("../../Data/TCITY_2015_EW_BGG_V2.shp") %>%
+# Town boundaries
+towns <- st_read("../Data/TCITY_2015_EW_BGG_V2.shp") %>%
   st_transform(4326) 
 towns$geometry <- st_make_valid(towns$geometry)
-towns_within_LCR_buffer <- towns[st_within(towns, LCR_bound_small_buffer, sparse = FALSE), ]
+towns_within_WYCA_buffer <- towns[st_within(towns, WYCA_bound_small_buffer, sparse = FALSE), ]
 
-# # Plot job locations levels across LCR -------------------------------------------
+## Plot job locations levels across WYCA -------------------------------------------
 # create job density variable
-LCR_dataset_jobs$job_density <- LCR_dataset_jobs$Employed_Population/LCR_dataset_jobs$LSOA_Area
+WYCA_dataset_jobs$job_density <- WYCA_dataset_jobs$Employed_Population/WYCA_dataset_jobs$LSOA_Area
 
 # Calculate PTJA quantile breaks
-breaks_jobs <- classIntervals(LCR_dataset_jobs$job_density, n = 5, style = "quantile")$brks
+breaks_jobs <- classIntervals(WYCA_dataset_jobs$job_density, n = 5, style = "quantile")$brks
 
 # Create a factor variable for PTJA coloring
-LCR_dataset_jobs$job_density_colour <- cut(LCR_dataset_jobs$job_density, breaks = breaks_jobs, include.lowest = TRUE, labels = FALSE)
-LCR_dataset_jobs$job_density_colour <- as.factor(LCR_dataset_jobs$job_density_colour)
-labels_jobs <- c("Lowest quintile","","","","Highest quintile")
-(jobs <- LCR_dataset_jobs %>%
+WYCA_dataset_jobs$job_density_colour <- cut(WYCA_dataset_jobs$job_density, 
+                                            breaks = breaks_jobs, 
+                                            include.lowest = TRUE, 
+                                            labels = FALSE)
+WYCA_dataset_jobs$job_density_colour <- as.factor(WYCA_dataset_jobs$job_density_colour)
+
+labels_jobs <- c("<150 jobs/km2","150 - 350","350 - 725","725 - 1650",">1650 jobs/km2")
+
+(jobs <- WYCA_dataset_jobs %>%
   ggplot() +
-  geom_sf(data=LCR_boundary, colour="black",linewidth=1.5)+
+  geom_sf(data=WYCA_boundary, colour="black",linewidth=1.5)+
   geom_sf(aes(fill = job_density_colour), col = NA) +
   scale_fill_brewer(palette = "Spectral", 
                     direction = -1,
                     labels = labels_jobs) +
-  labs(fill = "Job density (quintiles)")+
-    geom_sf(data=towns_within_LCR_buffer, fill = NA, col = "black", linewidth = 0.3) +
-    geom_sf(data=(towns_centres_LCR %>% filter("Buffer" == "N")), shape = 21, fill = 'white', size = 1.5) +
-    geom_label( x=-2.88, y=53.40, label="Liverpool", size=3) +
-    geom_label( x=-2.73, y=53.43, label="St. Helens", size=3) +
-    geom_label( x=-3.10, y=53.64, label="Southport", size=3) +
-    geom_label( x=-3.12, y=53.39, label="Birkenhead", size=3) +
+  labs(fill = "Job Density (quintiles)",
+       title = "West Yorkshire Combined Authority Job Distribution") +
+    theme(plot.title = element_text(hjust = 0.5)) +
+    geom_sf(data=towns_within_WYCA_buffer, fill = NA, col = "black", linewidth = 0.3) +
+    geom_sf(data=(towns_centres_WYCA %>% filter("Buffer" == "N")), shape = 21, fill = 'white', size = 1.5) +
+    geom_label( x=-1.756, y=53.79, label="Bradford", size=3) +
+    geom_label( x=-1.859, y=53.723, label="Halifax", size=3) +
+    geom_label( x=-1.782, y=53.646, label="Huddersfield", size=3) +
+    geom_label( x=-1.545, y=53.796, label="Leeds", size=3) +
+    geom_label( x=-1.5058, y=53.676, label="Wakefield", size=3)+
     annotation_north_arrow(  height = unit(1, "cm"),
                              width = unit(1, "cm")) +
     annotation_scale(location = "bl",
@@ -71,21 +77,4 @@ labels_jobs <- c("Lowest quintile","","","","Highest quintile")
                      pad_y = unit(0, "cm"))+
   theme_void() )
 
-ggsave(file = "Images/jobs_LCR.jpeg", device = "jpeg", plot = jobs)
-
-# Create job locations variable as % of LCR jobs
-LCR_dataset_jobs$Employed_Population_percent <- LCR_dataset_jobs$Employed_Population*100/sum(LCR_dataset_jobs$Employed_Population)
-# plot as % of total jobs in LCR
-LCR_dataset_jobs %>%
-  filter(job_density<50000)%>%
-  ggplot() +
-  geom_sf(data=LCR_boundary, colour="black",linewidth=1.5)+
-  geom_sf(aes(fill = job_density), col = NA) +
-  scale_fill_gradient(low="blue", high="green") +
-  labs(fill = "Job density
-       (jobs / sq.km)")+
-  annotation_north_arrow()+
-  annotation_scale(location = "bl",
-                   pad_x = unit(2.5, "cm"),
-                   pad_y = unit(0, "cm"))+
-  theme_void()
+ggsave(file = "WYCA/Images/jobs_WYCA.jpeg", device = "jpeg", plot = jobs)
diff --git a/WYCA/3.b.Population_Density_Map_WYCA.R b/WYCA/3.b.Population_Density_Map_WYCA.R
index 3b33484..be38528 100644
--- a/WYCA/3.b.Population_Density_Map_WYCA.R
+++ b/WYCA/3.b.Population_Density_Map_WYCA.R
@@ -1,4 +1,4 @@
-### RUN 3.c.PopulationDensityChange_Liverpool.R FIRST to create the accurate Population shapefile
+### RUN 3.c.PopulationDensityChange_WYCA.R FIRST to create the accurate Population shapefile
 # 1. Setup ----------------------------------------------------------------
 options(rgl.useNULL = FALSE)
 require(tidyverse)
@@ -15,11 +15,9 @@ require(magick)
 require(extrafont)
 library(stars)
 
-setwd("~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/Liverpool_City_Region")
-
 # 2. Load Data ------------------------------------------------------------
 # Read population density shapefile from "Population_Density_Change.R" script
-LCR_pop <- read_sf("../../Data/LCR_population.shp") %>%
+WYCA_pop <- read_sf("../Data/WYCA_population.shp") %>%
   st_transform(4326) %>%
   rename("Pop_Dens_change" = "Pp_dns_",
          "LSOA21CD" = "LSOA21C",
@@ -27,7 +25,7 @@ LCR_pop <- read_sf("../../Data/LCR_population.shp") %>%
   dplyr::select(LSOA21CD, Pop_Dens_change, Pop_Dens, geometry)
 
 # check the boundary plot
-ggplot(LCR_pop) +
+ggplot(WYCA_pop) +
   geom_sf(aes(fill = Pop_Dens),
           color = "gray66",
           linewidth = 0)+
@@ -35,9 +33,8 @@ ggplot(LCR_pop) +
 
 
 # Create Bounding Box -----------------------------------------------------
-
 # setting the boundary as a bounding box
-bbox <- st_bbox(LCR_pop) %>%
+bbox <- st_bbox(WYCA_pop) %>%
   st_set_crs(4326)
 
 # finding the aspect ratio
@@ -67,7 +64,7 @@ if(width > height) {
 size = 250 * 3.5
 
 pop_raster <- st_rasterize(
-  LCR_pop,
+  WYCA_pop,
   nx = floor(size * w_ratio) %>% as.numeric(),
   ny = floor(size * h_ratio) %>% as.numeric()
 )
@@ -127,7 +124,7 @@ render_camera(theta = 0,
 )
 
 
-outfile <- glue::glue("Plots/LCR_Pop.png")
+outfile <- glue::glue("WYCA/Images/WYCA_Pop.png")
 
 {
   start_time <- Sys.time()
diff --git a/WYCA/3.c.PopulationDensityChange_WYCA.R b/WYCA/3.c.PopulationDensityChange_WYCA.R
index ddd1980..6212f57 100644
--- a/WYCA/3.c.PopulationDensityChange_WYCA.R
+++ b/WYCA/3.c.PopulationDensityChange_WYCA.R
@@ -1,32 +1,29 @@
-
+## Script for mapping population density change in WYCA between 2001 and 2021
 # 1. Setup ----------------------------------------------------------------
 library(tidyverse)
 library(knitr)
 library(sf)
 
-setwd("~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/WYCA")
-
 # 2. Load and Filter Datasets ------------------------------------------------------------
-
 # WYCA Boundary + buffer
-Boundaries <- read_sf("../../Data/CAUTH_DEC_2023_EN_BFC.shp")
-WYCA_boundary <- Boundaries %>% filter(CAUTH23NM == "Liverpool City Region") %>%
+Boundaries <- read_sf("../Data/CAUTH_DEC_2023_EN_BFC.shp")
+WYCA_boundary <- Boundaries %>% filter(CAUTH23NM == "West Yorkshire") %>%
   st_transform(4326) 
 WYCA_bound_small_buffer <- WYCA_boundary %>% st_buffer(dist=1000)
 
 # Read Local Authority District (LAD) boundaries
-LADs <- read_sf("../../Data/LAD_DEC_2021_GB_BFC.shp") %>%
+LADs <- read_sf("../Data/LAD_DEC_2021_GB_BFC.shp") %>%
   st_transform(4326)
 
 # Town centres
-towns_centres_WYCA <- read_sf("../../Data/towns_centres_WYCA.shp") %>% 
+towns_centres_WYCA <- read_sf("../Data/towns_centres_WYCA.shp") %>% 
   st_transform(4326) %>%
   filter(buffer == "N") 
 
 # Filter LADs within WYCA
 LADs_WYCA <- LADs %>% filter(as.vector(st_within(., WYCA_bound_small_buffer, sparse = FALSE))) %>% 
   st_transform(4326)
-Pop_2001 <- read_csv("../../Data/Census2001_UsualPopulation.csv", skip=6) %>%
+Pop_2001 <- read_csv("../Data/Census2001_UsualPopulation.csv", skip=6) %>%
   rename("LSOA01" = "2001 super output areas - lower layer",
          "Pop_2001" = "2001") %>%
   drop_na() %>%
@@ -34,15 +31,23 @@ Pop_2001 <- read_csv("../../Data/Census2001_UsualPopulation.csv", skip=6) %>%
   mutate(Pop_2001 = as.numeric(Pop_2001))
 
 # Read Towns and City boundaries
-towns <- st_read("../../Data/TCITY_2015_EW_BGG_V2.shp") %>%
+towns <- st_read("../Data/TCITY_2015_EW_BGG_V2.shp") %>%
   st_transform(4326) 
 towns$geometry <- st_make_valid(towns$geometry)
 towns_within_WYCA <- towns[st_within(towns, WYCA_bound_small_buffer, sparse = FALSE), ]
 
+# Load 2021 Census data on population
+Pop2021 <- read_csv("../Data/Census2021_UsualPopulation.csv", skip=6) %>%
+  rename("LSOA21CD" = "2021 super output area - lower layer",
+         "Usual_Pop_2021" = "2021") %>%
+  drop_na() %>%
+  mutate(LSOA21CD = str_sub(LSOA21CD,1, 9))%>%
+  mutate(Usual_Pop_2021 = as.numeric(Usual_Pop_2021))
+
 
 # 2001 - 2011 LSOA Conversion --------------------------------------------------
 # Load 2001-2011 LSOA lookup table
-LSOA_lookup_01_11 <- read_csv("../../Data/LSOA_2001_2011_lookup.csv")
+LSOA_lookup_01_11 <- read_csv("../Data/LSOA_2001_2011_lookup.csv")
 
 # Split LSOAs - population split between new LSOAs generated in 2011
 # Join the LSOA lookup table to the dataset
@@ -79,9 +84,8 @@ Pop_2001_corrected_11 <- Pop_2001_corrected_11 %>%
 
 
 # 2011 - 2021 LSOA Conversion --------------------------------------------------
-
 # Load 2011-2021 LSOA lookup table
-LSOA_lookup_11_21 <- read_csv("../../Data/LSOA_(2011)_to_LSOA_(2021)_to_Local_Authority_District_(2022)_Lookup_for_England_and_Wales.csv")
+LSOA_lookup_11_21 <- read_csv("../Data/LSOA_(2011)_to_LSOA_(2021)_to_Local_Authority_District_(2022)_Lookup_for_England_and_Wales.csv")
 
 # Split LSOAs - population split between new LSOAs generated in 2021
 # Join the LSOA lookup table to the dataset
@@ -115,26 +119,24 @@ Pop_2001_corrected_21 <- Pop_2001_corrected_21 %>%
   left_join(summarized_population11, by = "LSOA21CD")  %>%
   mutate(Pop_2001_corrected_21 = coalesce(Pop_2001_corrected_11.y, Pop_2001_corrected_11.x)) %>%
   dplyr::select(c(Pop_2001_corrected_21, LSOA21CD))  
+# Step 4: replace NAs in Pop_2001_corrected_21 with population as stated in 2021. 
+# There are only 5 across England and Wales, but this prevents any NAs in the dataset
+Pop_2001_corrected_21 <- Pop2021 %>% 
+  left_join(Pop_2001_corrected_21, by = "LSOA21CD") %>%
+  mutate(Pop_2001_corrected_21 = coalesce(Pop_2001_corrected_21, Usual_Pop_2021)) %>%
+  dplyr::select(LSOA21CD, Pop_2001_corrected_21)
 
-
+# Calculate the scaling factor
+sum(Pop2021$Usual_Pop_2021, na.rm = TRUE)/sum(Pop_2001_corrected_21$Pop_2001_corrected_21, na.rm = TRUE)
 
 # Load West Yorkshire Dataset Shapefile ---------------------------------------
-
-WYCA_dataset <- read_sf("../../Data/WYCA_dataset.shp") %>%
+WYCA_dataset <- read_sf("../Data/WYCA_dataset.shp") %>%
                rename("LSOA21CD" = "LSOA21C",
                       "LSOA21NM" = "LSOA21N",
                       "PT_Job_Access_Index" = "PT_Jb_A_I",
                       "Employed_Population" = "Employd",
                       "LSOA_Area" = "LSOA__2") 
 
-# Load 2021 Census data on population
-Pop2021 <- read_csv("../../Data/Census2021_UsualPopulation.csv", skip=6) %>%
-  rename("LSOA21CD" = "2021 super output area - lower layer",
-         "Usual_Pop_2021" = "2021") %>%
-  drop_na() %>%
-  mutate(LSOA21CD = str_sub(LSOA21CD,1, 9))%>%
-  mutate(Usual_Pop_2021 = as.numeric(Usual_Pop_2021))
-sum(Pop2021$Usual_Pop_2021, na.rm = TRUE)/sum(Pop_2001_corrected_21$Pop_2001_corrected_21, na.rm = TRUE)
 
 # Join the population datasets to the WYCA dataset
 WYCA_population <- WYCA_dataset %>%
@@ -167,42 +169,43 @@ WYCA_population <- WYCA_population  %>%
 WYCA_population  %>%
       filter(as.vector(st_within(., WYCA_bound_small_buffer, sparse = FALSE))) %>%
           ggplot() +
-  #          geom_sf(data=WYCA_boundary, colour="black",linewidth=1.5) +
-            geom_sf(aes(fill = Pop_dens_change_cat), color = NA) +
+            geom_sf(data=WYCA_boundary, colour="black",linewidth=1.5) +
+            geom_sf(aes(fill = Pop_dens_change_cat), colour = NA) +
             scale_fill_brewer() +
             labs(fill=" Population Change per
             square kilometer",
                  title = "Population Density Change across WYCA from 2001-2021") +
- #           geom_sf(data=LADs_WYCA, fill = NA, col = "red", size = 1)  +
             geom_sf(data=towns_within_WYCA, fill = NA, col = "black", size = 1) +
             geom_sf(data=towns_centres_WYCA, shape = 21, fill = 'white', size = 1.5) +
-            geom_label( x=-2.88, y=53.40, label="Liverpool", size=3) +
-            geom_label( x=-2.73, y=53.43, label="St. Helens", size=3) +
-            geom_label( x=-3.10, y=53.64, label="Southport", size=3) +
-            geom_label( x=-3.12, y=53.39, label="Birkenhead", size=3) +
+            geom_label( x=-1.756, y=53.79, label="Bradford", size=3) +
+            geom_label( x=-1.859, y=53.723, label="Halifax", size=3) +
+            geom_label( x=-1.782, y=53.646, label="Huddersfield", size=3) +
+            geom_label( x=-1.545, y=53.796, label="Leeds", size=3) +
+            geom_label( x=-1.5058, y=53.676, label="Wakefield", size=3)+
             theme_void() 
 
-ggsave(file = "Images/Population_Density_Change_WYCA.jpg", device = "jpeg")
+ggsave(file = "WYCA/Images/Population_Density_Change_WYCA.jpg", device = "jpeg")
 
 # Write as output shapefile
-write_sf(WYCA_population, "../../Data/WYCA_population.shp")
+write_sf(WYCA_population, "../Data/WYCA_population.shp")
 
 # Calculate Pop_density quantile breaks
 breaks_pop <- classIntervals(WYCA_population$Pop_2021_density, n = 5, style = "quantile")$brks
 # Create a factor variable for PTJA coloring
 WYCA_population$Pop_density_cat <- cut(WYCA_population$Pop_2021_density, breaks = breaks_pop, include.lowest = TRUE, labels = FALSE)
 WYCA_population$Pop_density_cat <- as.factor(WYCA_population$Pop_density_cat)
-labels_pop <- c("<2,100","2,100 - 3,800","3,800 - 5,100","5,100 - 6,500","> 6,500")
+labels_pop <- c("<1,200","1,200 - 2,600","2,600 - 4,000","4,000 - 5,700","> 5,700")
 WYCA_population  %>%
   filter(as.vector(st_within(., WYCA_bound_small_buffer, sparse = FALSE))) %>%
   ggplot() +
-  #          geom_sf(data=WYCA_boundary, colour="black",linewidth=1.5) +
-  geom_sf(aes(fill = Pop_density_cat), color = NA) +
-  scale_fill_brewer(palette = "Spectral", 
-                    direction = -1,
-                    labels = labels_pop) +
-  labs(fill=" Population per square kilometer
-       (quintiles)",
-       title = "Population Density across WYCA in 2021") +
-  theme_void()
-ggsave(file = "Images/Population_Density_WYCA.jpg", device = "jpeg")
+      geom_sf(data=WYCA_boundary, colour="black",linewidth=1.5) +
+      geom_sf(aes(fill = Pop_density_cat), color = NA) +
+      scale_fill_brewer(palette = "Spectral", 
+                        direction = -1,
+                        labels = labels_pop) +
+      labs(fill=" Population per 
+        square kilometer
+           (quintiles)",
+           title = "Population Density across WYCA in 2021") +
+      theme_void()
+ggsave(file = "WYCA/Images/Population_Density_WYCA.jpg", device = "jpeg")
diff --git a/WYCA/3.d.PTJA_WYCA.R b/WYCA/3.d.PTJA_WYCA.R
index 040ac14..a5d97c9 100644
--- a/WYCA/3.d.PTJA_WYCA.R
+++ b/WYCA/3.d.PTJA_WYCA.R
@@ -9,12 +9,9 @@ library(ggplot2)
 library(ggspatial)
 library(classInt)
 
-setwd("~/Library/CloudStorage/GoogleDrive-sam.allwood3@gmail.com/My Drive/Consulting/Unemployment_Public_Transport_Access/Liverpool_City_Region")
-
 # 2. Load Data ------------------------------------------------------------
-
-# Load Liverpool Geo data from file
-LCR_dataset_PTJA <- read_sf("../../Data/LCR_dataset.shp") %>%
+# Load WYCA Geo data from file
+WYCA_dataset_PTJA <- read_sf("../Data/WYCA_dataset.shp") %>%
   rename("LSOA_Code" = "LSOA21C",
          "LSOA_Name" = "LSOA21N",
          "PT_Job_Access_Index" = "PT_Jb_A_I",
@@ -30,56 +27,57 @@ LCR_dataset_PTJA <- read_sf("../../Data/LCR_dataset.shp") %>%
          PT_Job_Access_Index_train = PT_Job_Access_Index_train - PT_Job_Access_Index_walk) 
 
 
-# LCR Boundary + buffer
-Boundaries <- read_sf("../../Data/CAUTH_DEC_2023_EN_BFC.shp")
-LCR_boundary <- Boundaries %>% filter(CAUTH23NM == "Liverpool City Region") %>%
+# WYCA Boundary + buffer
+Boundaries <- read_sf("../Data/CAUTH_DEC_2023_EN_BFC.shp")
+WYCA_boundary <- Boundaries %>% filter(CAUTH23NM == "West Yorkshire") %>%
   st_transform(4326) 
-LCR_bound_small_buffer <- LCR_boundary %>% st_buffer(dist=25)
+WYCA_bound_small_buffer <- WYCA_boundary %>% st_buffer(dist=25)
 
 # Town centres
-towns_centres_LCR <- read_sf("../../Data/towns_centres_LCR.shp") %>%
+towns_centres_WYCA <- read_sf("../Data/towns_centres_WYCA.shp") %>%
   filter(buffer == "N")
 
-# # Plot PTJA index in LCR ---------------------------------------------
-LCR_dataset_PTJA %>%
+# # Plot PTJA index in WYCA ---------------------------------------------
+WYCA_dataset_PTJA %>%
   ggplot() +
-  geom_sf(data=LCR_boundary, colour="red",linewidth=1.5)+
+  geom_sf(data=WYCA_boundary, colour="red",linewidth=1.5)+
   geom_sf(aes(fill = PT_Job_Access_Index), col = NA) +
   scale_fill_gradient(low="BLACK", high="WHITE") +
   labs(fill = "PTJA Index")+
   theme_void()
 
 
-# Plot PTJA index in LCR by quantile
+# Plot PTJA index in WYCA by quantile
 # Calculate PTJA quantile breaks
-breaks_PTJA <- classIntervals(LCR_dataset_PTJA$PT_Job_Access_Index, n = 5, style = "quantile")$brks
+breaks_PTJA <- classIntervals(WYCA_dataset_PTJA$PT_Job_Access_Index, n = 5, style = "quantile")$brks
 breaks_PTJA[1] <- 00
 # Create a factor variable for PTJA coloring
-LCR_dataset_PTJA$color_PTJA <- cut(LCR_dataset_PTJA$PT_Job_Access_Index, breaks = breaks_PTJA, include.lowest = TRUE, labels = FALSE)
-LCR_dataset_PTJA$color_PTJA <- as.factor(LCR_dataset_PTJA$color_PTJA)
-LCR_dataset_PTJA$color_PTJA_bus <- cut(LCR_dataset_PTJA$PT_Job_Access_Index_bus, breaks = breaks_PTJA, include.lowest = TRUE, labels = FALSE)
-LCR_dataset_PTJA$color_PTJA_bus <- as.factor(LCR_dataset_PTJA$color_PTJA_bus)
-LCR_dataset_PTJA$color_PTJA_train <- cut(LCR_dataset_PTJA$PT_Job_Access_Index_train, breaks = breaks_PTJA, include.lowest = TRUE, labels = FALSE)
-LCR_dataset_PTJA$color_PTJA_train <- as.factor(LCR_dataset_PTJA$color_PTJA_train)
+WYCA_dataset_PTJA$color_PTJA <- cut(WYCA_dataset_PTJA$PT_Job_Access_Index, breaks = breaks_PTJA, include.lowest = TRUE, labels = FALSE)
+WYCA_dataset_PTJA$color_PTJA <- as.factor(WYCA_dataset_PTJA$color_PTJA)
+WYCA_dataset_PTJA$color_PTJA_bus <- cut(WYCA_dataset_PTJA$PT_Job_Access_Index_bus, breaks = breaks_PTJA, include.lowest = TRUE, labels = FALSE)
+WYCA_dataset_PTJA$color_PTJA_bus <- as.factor(WYCA_dataset_PTJA$color_PTJA_bus)
+WYCA_dataset_PTJA$color_PTJA_train <- cut(WYCA_dataset_PTJA$PT_Job_Access_Index_train, breaks = breaks_PTJA, include.lowest = TRUE, labels = FALSE)
+WYCA_dataset_PTJA$color_PTJA_train <- as.factor(WYCA_dataset_PTJA$color_PTJA_train)
 labels_PTJA <- c("Lowest PTJA Quintile","","","","Highest PTJA Quintile")
 
-# Plot PTJA by quintile in LCR
-(PTJA <- LCR_dataset_PTJA %>%
+# Plot PTJA by quintile in WYCA
+(PTJA <- WYCA_dataset_PTJA %>%
     dplyr::select(color_PTJA, geometry) %>%
   ggplot() +
-  geom_sf(data=LCR_boundary, colour="black",linewidth=1.0) +
-  geom_sf(data = LCR_dataset_PTJA, aes(fill = color_PTJA), color = NA) + # color = NA removes the LSOA boundaries
+  geom_sf(data=WYCA_boundary, colour="black",linewidth=1.0) +
+  geom_sf(data = WYCA_dataset_PTJA, aes(fill = color_PTJA), color = NA) + # color = NA removes the LSOA boundaries
   scale_fill_brewer(palette = "Spectral", 
                     direction = -1,
                     labels = labels_PTJA) +
-  labs(fill = "Public Transport
-       Job Accessibility 
-       (quintiles)") +
-  geom_sf(data=towns_centres_LCR, shape = 21, fill = 'white', size = 1.5) +
-    geom_label( x=-2.88, y=53.40, label="Liverpool", size=3) +
-    geom_label( x=-2.73, y=53.43, label="St. Helens", size=3) +
-    geom_label( x=-3.10, y=53.64, label="Southport", size=3) +
-    geom_label( x=-3.12, y=53.39, label="Birkenhead", size=3) +
+  labs(fill = "PTJA 
+       (quintiles)",
+       title = "Public Transport Job Accessibility in West Yorkshire Combined Authority") +
+  geom_sf(data=towns_centres_WYCA, shape = 21, fill = 'white', size = 1.5) +
+    geom_label( x=-1.756, y=53.83, label="Bradford", size=3) +
+    geom_label( x=-1.859, y=53.75, label="Halifax", size=3) +
+    geom_label( x=-1.782, y=53.675, label="Huddersfield", size=3) +
+    geom_label( x=-1.545, y=53.825, label="Leeds", size=3) +
+    geom_label( x=-1.5058, y=53.70, label="Wakefield", size=3)+
     annotation_north_arrow(height = unit(1, "cm"),
                            width = unit(1, "cm"),) +
     annotation_scale(location = "bl",
@@ -87,14 +85,14 @@ labels_PTJA <- c("Lowest PTJA Quintile","","","","Highest PTJA Quintile")
                      pad_y = unit(0, "cm"))+
   theme_void())
 
-ggsave(file = "Images/PTJA.jpeg", device = "jpeg", plot = PTJA)
+ggsave(file = "WYCA/Images/PTJA.jpeg", device = "jpeg", plot = PTJA)
 
 # bus contribution map
-(PTJA_bus <- LCR_dataset_PTJA %>%
+(PTJA_bus <- WYCA_dataset_PTJA %>%
     dplyr::select(color_PTJA_bus, geometry) %>%
     ggplot() +
-    geom_sf(data=LCR_boundary, colour="black",linewidth=1.0) +
-    geom_sf(data = LCR_dataset_PTJA, aes(fill = color_PTJA_bus), color = NA) + # color = NA removes the LSOA boundaries
+    geom_sf(data=WYCA_boundary, colour="black",linewidth=1.0) +
+    geom_sf(data = WYCA_dataset_PTJA, aes(fill = color_PTJA_bus), color = NA) + # color = NA removes the LSOA boundaries
     scale_fill_brewer(palette = "Spectral", 
                       direction = -1,
                       labels = c("Lowest",
@@ -104,15 +102,15 @@ ggsave(file = "Images/PTJA.jpeg", device = "jpeg", plot = PTJA)
                                  "Highest")) +
     labs(fill = "Public Transport Job 
          Accessibility: Bus")+
-    geom_sf(data=towns_centres_LCR, shape = 21, fill = 'white', size = 1.5) +
+    geom_sf(data=towns_centres_WYCA, shape = 21, fill = 'white', size = 1.5) +
     theme_void())
 
 # tram contribution map
-(PTJA_train <- LCR_dataset_PTJA %>%
+(PTJA_train <- WYCA_dataset_PTJA %>%
     dplyr::select(color_PTJA_train, geometry) %>%
     ggplot() +
-    geom_sf(data=LCR_boundary, colour="black",linewidth=1.0) +
-    geom_sf(data = LCR_dataset_PTJA, aes(fill = color_PTJA_train), color = NA) + # color = NA removes the LSOA boundaries
+    geom_sf(data=WYCA_boundary, colour="black",linewidth=1.0) +
+    geom_sf(data = WYCA_dataset_PTJA, aes(fill = color_PTJA_train), color = NA) + # color = NA removes the LSOA boundaries
     scale_fill_brewer(palette = "Spectral", 
                       direction = -1,
                       labels = c("Lowest","","","","Highest")) +
@@ -123,12 +121,12 @@ ggsave(file = "Images/PTJA.jpeg", device = "jpeg", plot = PTJA)
 
 # Adjusted for Demand Potential -------------------------------------------
 # Calculate PTJA-D quantile breaks
-breaks_PTJA_D <- classIntervals(LCR_dataset_PTJA$PT_Job_Access_Index_Demand, 
+breaks_PTJA_D <- classIntervals(WYCA_dataset_PTJA$PT_Job_Access_Index_Demand, 
                                 n = 5, 
                                 style = "quantile")$brks
 breaks_PTJA_D[1] <- 00
 # Create a factor variable for PTJA coloring
-LCR_dataset_PTJA$color_PTJA_D <- cut(LCR_dataset_PTJA$PT_Job_Access_Index_Demand, 
+WYCA_dataset_PTJA$color_PTJA_D <- cut(WYCA_dataset_PTJA$PT_Job_Access_Index_Demand, 
                                        breaks = breaks_PTJA_D, 
                                        include.lowest = TRUE, 
                                        labels = FALSE) %>%
@@ -137,24 +135,25 @@ LCR_dataset_PTJA$color_PTJA_D <- cut(LCR_dataset_PTJA$PT_Job_Access_Index_Demand
 
 labels_PTJA_D <- c("Lowest PTJA-D Quintile","","","","Highest PTJA-D Quintile")
 
-# Plot PTJA by quintile in LCR
-(PTJA_D <- LCR_dataset_PTJA %>%
+# Plot PTJA by quintile in WYCA
+(PTJA_D <- WYCA_dataset_PTJA %>%
     dplyr::select(color_PTJA_D, geometry) %>%
     ggplot() +
-    geom_sf(data=LCR_boundary, colour="black",linewidth=1.0) +
-    geom_sf(data = LCR_dataset_PTJA, aes(fill = color_PTJA_D), color = NA) + # color = NA removes the LSOA boundaries
+    geom_sf(data=WYCA_boundary, colour="black",linewidth=1.0) +
+    geom_sf(data = WYCA_dataset_PTJA, aes(fill = color_PTJA_D), color = NA) + # color = NA removes the LSOA boundaries
     scale_fill_brewer(palette = "Spectral", 
                       direction = -1,
                       labels = labels_PTJA_D) +
-    labs(fill = "Public Transport
-       Job Accessibility 
-       (Demand-adjusted)
-       (quintiles)")+
-    geom_sf(data=towns_centres_LCR, shape = 21, fill = 'white', size = 1.5) +
-    geom_label( x=-2.88, y=53.40, label="Liverpool", size=3) +
-    geom_label( x=-2.73, y=53.43, label="St. Helens", size=3) +
-    geom_label( x=-3.10, y=53.64, label="Southport", size=3) +
-    geom_label( x=-3.12, y=53.39, label="Birkenhead", size=3) +
+    labs(fill = "PTJA-D
+       (quintiles)",
+         title = "Public Transport Job Accessibility (Demand-adjusted) in 
+         West Yorkshire Combined Authority") +
+    geom_sf(data=towns_centres_WYCA, shape = 21, fill = 'white', size = 1.5) +
+    geom_label( x=-1.756, y=53.83, label="Bradford", size=3) +
+    geom_label( x=-1.859, y=53.75, label="Halifax", size=3) +
+    geom_label( x=-1.782, y=53.675, label="Huddersfield", size=3) +
+    geom_label( x=-1.545, y=53.825, label="Leeds", size=3) +
+    geom_label( x=-1.5058, y=53.70, label="Wakefield", size=3)+
     annotation_north_arrow( height = unit(1, "cm"),
                             width = unit(1, "cm")) +
     annotation_scale(location = "bl",
@@ -162,4 +161,4 @@ labels_PTJA_D <- c("Lowest PTJA-D Quintile","","","","Highest PTJA-D Quintile")
                      pad_y = unit(0, "cm")) +
     theme_void())
 
-ggsave(file = "Images/PTJA_D.jpeg", device = "jpeg", plot = PTJA_D)
+ggsave(file = "WYCA/Images/PTJA_D.jpeg", device = "jpeg", plot = PTJA_D)
diff --git a/WYCA/Images/PTJA.jpeg b/WYCA/Images/PTJA.jpeg
new file mode 100644
index 0000000..4222a31
Binary files /dev/null and b/WYCA/Images/PTJA.jpeg differ
diff --git a/WYCA/Images/PTJA_D.jpeg b/WYCA/Images/PTJA_D.jpeg
new file mode 100644
index 0000000..3452694
Binary files /dev/null and b/WYCA/Images/PTJA_D.jpeg differ
diff --git a/WYCA/Images/Population_Density_Change_WYCA.jpg b/WYCA/Images/Population_Density_Change_WYCA.jpg
new file mode 100644
index 0000000..b3e4baf
Binary files /dev/null and b/WYCA/Images/Population_Density_Change_WYCA.jpg differ
diff --git a/WYCA/Images/Population_Density_WYCA.jpg b/WYCA/Images/Population_Density_WYCA.jpg
new file mode 100644
index 0000000..4ca66a0
Binary files /dev/null and b/WYCA/Images/Population_Density_WYCA.jpg differ
diff --git a/WYCA/Images/WYCA_Pop.png b/WYCA/Images/WYCA_Pop.png
new file mode 100644
index 0000000..74bc433
Binary files /dev/null and b/WYCA/Images/WYCA_Pop.png differ
diff --git a/WYCA/Images/jobs_WYCA.jpeg b/WYCA/Images/jobs_WYCA.jpeg
new file mode 100644
index 0000000..1eaa4ee
Binary files /dev/null and b/WYCA/Images/jobs_WYCA.jpeg differ