#30DayMapChallenge 2021 - Day 07 - Green

street_grid.R

# Street orientation grids of New Zealand cities
# for #30DayMapChallenge 2021 - Day 07 - Green
#  -- David Friggens, November 2021

# Inspired by:
# Geoff Boeing
#  https://geoffboeing.com/2018/07/comparing-city-street-orientations/
#  https://geoffboeing.com/2019/09/urban-street-network-orientation/
# Erin Davis 
#  https://erdavis.com/2019/11/21/street-network-orientation-by-road-type/


# Setup -------------------------------------------------------------------

library(magrittr)
library(readr)
library(dplyr)
library(purrr)
library(sf)
library(stringr)
library(geosphere)
library(ggplot2)
library(patchwork)

PLOT_FONT <- "Nunito Sans"
# extracted from images on nzta.govt.nz
SIGN_GREEN <- "#00973f"
BG_COLOUR <- "white"


# Read data --------------------------------------------------------------

# Table 7981 from https://nzdotstat.stats.govt.nz
population <- 
  read_tsv("population/population_ur.tsv") %>% 
  filter(YEAR == 2021) %>% 
  transmute(
    area_code = AREA %>% str_replace("^URA", ""),
    area_name = Area,
    population = `Value  Flags`
  ) %>% 
  arrange(desc(population))

# https://data.linz.govt.nz/layer/53382-nz-roads-addressing/
roads <- 
  read_sf("roads/nz-roads-addressing.gpkg")

# https://datafinder.stats.govt.nz
urban_areas <- 
  read_sf("statsnz/urban-rural-2021-clipped-generalised.gpkg") %>% 
  select(area_code = UR2021_V1_00,
         area_name = UR2021_V1_00_NAME,
         urban_type = IUR2021_V1_00_NAME) %>% 
  filter(urban_type %>% str_detect("rban"))

roads <- 
  roads %>% 
  st_transform(crs = st_crs(urban_areas))



# Define functions -------------------------------------------------------------------

get_first_coordinates <- 
  function(my_linestring) {
    my_linestring %>% 
      st_cast("POINT") %>% 
      head(1) %>% 
      st_transform(crs = 4193L) %>% 
      st_coordinates()
  }

get_last_coordinates <- 
  function(my_linestring) {
    my_linestring %>% 
      st_cast("POINT") %>% 
      tail(1) %>% 
      st_transform(crs = 4193L) %>% 
      st_coordinates()
  }



plot_my_roads <- 
  function(my_urban_code) {
    my_urban_name <- 
      urban_areas %>% 
      st_drop_geometry() %>% 
      filter(area_code == my_urban_code) %>% 
      pull(area_name)
    my_urban_roads <- 
      urban_areas %>% 
      filter(area_code == my_urban_code) %>% 
      st_intersection(roads) %>% 
      select(geom) %>% 
      # It's currently a mix, so going direct to LS only takes the first LT from
      # every MLS row. But if everything is MLS then they're expanded to create
      # rows for every LS
      st_cast("MULTILINESTRING") %>% 
      st_cast("LINESTRING")
    
    my_urban_grid <- 
      my_urban_roads %>% 
      rowwise() %>% 
      mutate(road_start = geom %>% get_first_coordinates(),
             road_end = geom %>% get_last_coordinates()) %>% 
      st_drop_geometry() %>% 
      mutate(road_bearing =
               bearing(road_start, road_end),
             road_distance =
               distHaversine(road_start, road_end)) %>% 
      select(road_bearing, road_distance) %>% 
      ungroup() %>% 
      mutate(road_bearing = road_bearing %>% round(-1))
    my_urban_grid <- 
      bind_rows(
        my_urban_grid,
        my_urban_grid %>% 
          mutate(road_bearing =
                   if_else(road_bearing >= 180,
                           road_bearing - 180,
                           road_bearing + 180))
      ) %>% 
      mutate(road_bearing =
               if_else(road_bearing < 0,
                       road_bearing + 360,
                       road_bearing) %>%
               subtract(5))
    
    my_urban_plot <- 
      ggplot(my_urban_grid,
           aes(x = road_bearing,
               weight = road_distance)) +
      geom_histogram(binwidth = 10, boundary = 5,
                     size = .01, closed = "left",
                     fill = SIGN_GREEN,
                     colour = "black") +
      scale_x_continuous(breaks = c(0, 90, 180, 270),
                         limits = c(-5, 355),
                         labels = c("N", "E", "S", "W")) +
      coord_polar(start = -pi/36) +
      xlab(NULL) + ylab(NULL) +
      ggtitle(my_urban_name) +
      theme_bw(base_family = PLOT_FONT) +
      theme(plot.title = element_text(family = PLOT_FONT, hjust = 0.5,
                                      size = 32,
                                      colour = SIGN_GREEN),
            axis.title.y = element_blank(),
            axis.text.y = element_blank(),
            axis.ticks.y = element_blank(),
            axis.text.x = element_text(colour = SIGN_GREEN),
            panel.border = element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank())
    
    return(my_urban_plot)
  }



# Run functions -----------------------------------------------------------

# So, this may not be super efficient - it took over an hour to run!
# But it came out fine and I don't have time to refactor now. :-)
urban_grids <- 
  population %>%
  head(20) %>% 
  mutate(grids = area_code %>% map(plot_my_roads))

gg_urban <- 
  urban_grids %>% 
  pull(grids) %>% 
  wrap_plots(ncol = 5) +
  plot_annotation(title = "Street Orientations: New Zealand Edition",
                  subtitle = "#30DayMapChallenge 2021 - Day 7 - Green",
                  caption = "David Friggens @dakvid, Data: Toit\u16b Te Whenua LINZ",
                  theme = theme(plot.title = element_text(size = 72, hjust = 0.5,
                                                          family = PLOT_FONT, colour = SIGN_GREEN),
                                plot.subtitle = element_text(size = 32, hjust = 0.5,
                                                             family = PLOT_FONT, colour = SIGN_GREEN),
                                plot.caption = element_text(size = 32,
                                                            family = PLOT_FONT, colour = SIGN_GREEN))) &
  theme(plot.background = element_rect(fill = BG_COLOUR, colour = BG_COLOUR),
        panel.background = element_rect(fill = BG_COLOUR, colour = BG_COLOUR))

ggsave(plot = gg_urban,
       path = "Day_07", filename = "Day_07_street_orientation.png", device = "png",
       width = 30, height = 30, dpi = 72)