new_sites.r

# Aim: generate scenarios for new developments

library(tidyverse)
library(sf)
library(stplanr)

# set-up and parameters ---------------------------------------------------

# setwd("~/cyipt/actdev") # run this script from the actdev folder
if(!exists("site_name")) { # assume all presets loaded if site_name exists
  site_name = "great-kneighton"   # which site to look at (can change)
  data_dir = "data-small" # for test sites
  max_length = 20000 # maximum length of desire lines in m
  household_size = 2.3 # mean UK household size at 2011 census
  min_flow_routes = 10 # threshold above which OD pairs are included
  region_buffer_dist = 2000
  large_area_buffer = 500
}

if(!exists("centroids_msoa")) {
  # run the build script if national data is missing
  source("code/build-setup.R")
}

# to load a custom site
site = sf::read_sf("new-site-demo/exeter-red-cow-village.geojson")

site_name = "exeter-red-cow-village"
site$site_name = site_name
site$full_name = "Exeter Red Cow Village (Liveable Exeter)"
site$dwellings_when_complete = 664

msoa_pops = readxl::read_xls(path = "data/mid2011msoaunformattedfile.xls", sheet = "Mid-2011 Persons", )
msoa_pops = msoa_pops %>% 
  select(geo_code1 = Code, msoa_population = "All Ages")

# estimated site populations
site_pops = sites %>% 
  st_drop_geometry() %>% 
  mutate(site_population = dwellings_when_complete * household_size)

message("Building for ", site$site_name)

path = file.path(data_dir, site_name)

zones_touching_site = zones_msoa_national[site, , op = sf::st_intersects]

zones_touching_site$overlap_size = units::drop_units(st_area(zones_touching_site))
zones_touching_site = zones_touching_site %>% 
  filter(overlap_size > 10000) %>% 
  select(-overlap_size)

# Route from site centroid (rather than MSOA centroid) --------------------
# `disaggregate.R` changes this to route from a random selection of homes within the site, to better represent the accessibility of the site as a whole
site_centroid = site %>% 
  st_transform(27700) %>% 
  st_centroid() %>% 
  st_transform(4326)

zone_data = zones_touching_site %>% 
  st_drop_geometry() %>%
  mutate(site_name = site$site_name)
site_c = right_join(site_centroid, zone_data) %>%
  select(geo_code, site_name)

# Generate desire lines ---------------------------------------------------
od_site = od %>% 
  filter(geo_code1 %in% zones_touching_site$geo_code) %>% 
  filter(geo_code2 %in% centroids_msoa$msoa11cd)
# intra-zonal flows are included, with the desire line going from the site centroid to the msoa centroid
# intra-zonal flows could later be represented using more detailed od-workplace zone data.

desire_lines_site = od::od_to_sf(x = od_site, z = site_c, zd = centroids_msoa)
desire_lines_site = desire_lines_site %>% 
  mutate(site_name = site_name)

# Adjust flows to represent site population, not MSOA population(s) -------
# for both MSOAs and development sites, these are entire populations, not commuter populations
desire_lines_site = inner_join(desire_lines_site, msoa_pops)
# desire_lines_site = inner_join(desire_lines_site, site_pops)
desire_lines_site = desire_lines_site %>% 
  mutate(site_population = site$dwellings_when_complete * 2.4)
site_population = unique(desire_lines_site$site_population)

unique_msoa_pops = desire_lines_site %>% 
  st_drop_geometry() %>% 
  select(geo_code1, msoa_population) %>%
  unique()
sum_msoa_pops = sum(unique_msoa_pops$msoa_population)
desire_lines_site = desire_lines_site %>% 
  mutate(sum_msoa_pops = sum_msoa_pops)

# keeping converted flows in the original columns
desire_lines_pops = desire_lines_site %>% 
  mutate(across(all:other, .fns = ~ ./ sum_msoa_pops * site_population))

# todo: add empirical data on 'new homes' effect (residents of new homes are more likely to drive than residents of older homes)
# could also adjust the base walking and cycling mode shares in response to the difference between journey distance from the site centroid as compared to journey distance from the MSOA centroid (eg in Cambridge, the MSOA centroid is a fair bit closer to the city centre than the site centroid, which could explain why such a high proportion of commuters are shown walking to work in the city centre)  

# For sites with 2 or more origin MSOAs, combine flows to avoid having multiple desire lines to the same destination MSOA
desire_lines_combined = desire_lines_pops %>% 
  group_by(geo_code2) %>% 
  summarise(
    geo_code1 = geo_code1[1], # do we even need this?
    across(all:other, sum)
  )

desire_lines_combined$length = round(stplanr::geo_length(desire_lines_combined))

# todo: add PT
desire_lines_combined = desire_lines_combined %>% 
  mutate(
    trimode_base = foot + bicycle + car_driver,
    pwalk_base = foot/trimode_base,
    pcycle_base = bicycle/trimode_base,
    pdrive_base = car_driver/trimode_base
  )
desire_lines_combined[is.na(desire_lines_combined)] = 0

desire_lines_combined = desire_lines_combined %>% 
  select(geo_code1, geo_code2, all, trimode_base, from_home:other, length, pwalk_base:pdrive_base) %>%
  mutate(across(where(is.numeric), round, 6))

st_precision(desire_lines_combined) = 1000000

dsn = file.path(data_dir, site_name, "all-census-od.csv")
obj = st_drop_geometry(desire_lines_combined)
if(file.exists(dsn)) file.remove(dsn)
if(!dir.exists(path)) dir.create(path)
readr::write_csv(obj, file = dsn)

# Round decimals and select sets of desire lines --------------------------
desire_lines_rounded = desire_lines_combined %>% 
  rename(all_base = all, walk_base = foot, cycle_base = bicycle, drive_base = car_driver) %>%
  mutate(
    across(all_base:other, smart.round),
    trimode_base = walk_base + cycle_base + drive_base
  ) %>% 
  filter(trimode_base > 0)

desire_lines_20km = desire_lines_rounded %>% 
  filter(length <= max_length)

desire_lines_threshold = desire_lines_rounded %>%
  filter(trimode_base >= min_flow_routes)

desire_lines_bounding = desire_lines_20km %>% 
  filter(geo_code2 %in% desire_lines_threshold$geo_code2)

# Large study area MSOAs --------------------------------------------------
large_study_area = sf::st_convex_hull(sf::st_union(desire_lines_bounding))
large_study_area = stplanr::geo_buffer(large_study_area, dist = large_area_buffer)

dsn = file.path(data_dir, site_name, "large-study-area.geojson")
if(file.exists(dsn)) file.remove(dsn)
sf::write_sf(large_study_area, dsn = dsn)

desire_lines_many = desire_lines_rounded[large_study_area, , op = sf::st_within]

desire_lines_many = desire_lines_many %>% 
  select(geo_code1, geo_code2, all_base, trimode_base, walk_base, cycle_base, drive_base, length, pwalk_base:pdrive_base)

# is this meant to be wrapped in {}?
if(!exists("disaggregate_desire_lines"))
  disagreggate_desire_lines = FALSE

if(disaggregate_desire_lines) {
  zones_many = zones_msoa_national %>% filter(geo_code %in% desire_lines_many$geo_code2)
  centroids_lsoa_national = pct::get_pct(layer = "c", national = TRUE)
  zones_lsoa_national = pct::get_pct(layer = "z", national = TRUE)
  centroids_lsoa_many = centroids_lsoa_national[zones_many, ]
  zones_lsoa_many = centroids_lsoa_national %>% filter(geo_code %in% centroids_lsoa_many$geo_code)
  desire_lines_many_min = desire_lines_many %>% select(geo_code1:drive_base) %>% sf::st_drop_geometry()
  n_lines = max(nrow(desire_lines_many), 20) # minimum of 20 desire lines
  p = sum(desire_lines_many_min$all_base) / n_lines
  # desire_lines_many_commute = od::od_disaggregate(od = desire_lines_many_min, z = zones_many, subzones = zones_lsoa_many) 
  file.edit("~/itsleeds/od/R/aggregate.R") # edit locally for testing + debugging
  desire_lines_disag = od_disaggregate(od = desire_lines_many_min, z = zones_many, subzones = zones_lsoa_many, population_per_od = p)
  
  # sanity tests  
  sum(desire_lines_many$all_base) == sum(desire_lines_disag$all_base)
  sum(desire_lines_many$walk_base) == sum(desire_lines_disag$walk_base)
  
  desire_lines_disag = desire_lines_disag %>% 
    mutate(
      trimode_base = foot + bicycle + car_driver,
      pwalk_base = foot/trimode_base,
      pcycle_base = bicycle/trimode_base,
      pdrive_base = car_driver/trimode_base
    )
  
}

# Create routes and generate Go Dutch scenario ---------------------
obj = desire_lines_many %>% select(-length) %>% top_n(n = 3, wt = all_base)
obj2 = desire_lines_many %>% filter(length < 6000) %>% select(-length) %>% top_n(n = 3, wt = all_base)

routes_fast = stplanr::route(l = obj, route_fun = cyclestreets::journey)
routes_balanced = stplanr::route(l = obj, route_fun = cyclestreets::journey, plan = "balanced")
routes_quiet = stplanr::route(l = obj, route_fun = cyclestreets::journey, plan = "quietest")