## The Red Rooster line

## Initiated by a radio interview described on https://www.abc.net.au/news/2023-11-25/disadvantaged-students-medicine-university/103141050

library(osmdata)
library(sf)
library(sfext) # Used to quickly calculate aspect ratio of bounding box. Install with pak::pkg_install("elipousson/sfext")
library(maptiles)
library(tidyterra)
library(ggplot2)
library(pals)
library(ggimage)
library(magick)

## Obtain boundaries for significant urban areas in Australia
## https://www.abs.gov.au/statistics/standards/australian-statistical-geography-standard-asgs-edition-3/jul2021-jun2026/access-and-downloads/digital-boundary-files

shapeAustralia <- read_sf("~/Downloads/SUA_2021_AUST_GDA94.shp")

## Create bounding boxes around the city shapes
bb <- list()
bb[["Perth"]] <- st_bbox(st_transform(st_as_sf(shapeAustralia[shapeAustralia$SUA_NAME21 == "Perth",]), 4326))
bb[["Sydney"]] <- st_bbox(st_as_sf(shapeAustralia[shapeAustralia$SUA_NAME21 == "Sydney",]))

## Obtain OpenStreetMap data on the location of fast food restaurants across the cities

fastfood <- lapply(bb, function(x) {
    opq(x, timeout = 50) |>
    add_osm_feature("amenity", "fast_food") |>
    add_osm_feature("access", "!private") |>
    osmdata_sf()})

## Map everything using downloaded background tiles

bbsfPerth <- st_as_sf(shapeAustralia[shapeAustralia$SUA_NAME21 == "Perth",])
 TODO - find B/W tile server
mapPerth <- get_tiles(x = bbsfPerth, provider = "OpenStreetMap", zoom = 9)
dev.new(width = 15, height = 15 * 9/16, unit = "cm")
ggplot(fastfood[["Perth"]]$osm_points) +
    geom_spatraster_rgb(data = mapPerth) +
    geom_sf(aes(colour = brand),
            size = 2) +
    theme_void() +
    scale_color_manual(values=as.vector(glasbey())) +
    theme(legend.position="right") +
    labs(title = "Fastfood in Perth metro",
         colour = "Brand")

## Obtain way data from OpenStreetMap in order to construct our own background map

majorRoads <- lapply(bb, function(x) {
    opq(x, timeout = 120) |>
    add_osm_feature(key = "highway", 
                  value = c("motorway", "primary", "secondary")) |>
    osmdata_sf()})

minorRoads <- lapply(bb, function(x) {
    opq(x, timeout = 120) |>
        add_osm_feature(key = "highway", value = c("tertiary")) |>
        osmdata_sf()})

boundaries <- lapply(bb, function(x) {
    opq(x, timeout = 120) |>
        add_osm_feature(key = "boundary", value = c("administrative")) |>
        add_osm_feature(key = "admin_level", value = c(2,8,9,10,11)) |>
        osmdata_sf()
})

water <- lapply(bb, function(x) {
    opq(x, timeout = 180) |>
    add_osm_feature(key = "natural", value = "water") |>
    add_osm_feature(key = 'name') |>
    osmdata_sf()})

## Create city outlines, sfc bounding box and separate out ocean area as OSM does not have that polygon

outlines <- sapply(names(bb), function(x) {
    st_union(boundaries[[x]]$osm_multipolygons[which(boundaries[[x]]$osm_multipolygons$admin_level == 9), ])},
    simplify = FALSE, USE.NAMES = TRUE)
bbSfc <- sapply(names(bb), function(x) {
    a <- st_as_sfc(bb[[x]])
    a <- st_transform(a, crs = st_crs(outlines[[x]]))
    return(a)},
    simplify = FALSE, USE.NAMES = TRUE)
ocean <- sapply(names(bbSfc), function(x) {
    a <- st_difference(bbSfc[[x]],outlines[[x]])
    a <- st_transform(a, crs = st_crs(outlines[[x]]))
    return(a)},
    simplify = FALSE, USE.NAMES = TRUE)

## Chop all the osm data to the sfc bounding box, as the queried data runs over its bounding box for many features

boundariesChopped <- sapply(names(bbSfc), function(x) {
    st_intersection(boundaries[[x]]$osm_multipolygons[which(boundaries[[x]]$osm_multipolygons$admin_level == 9), ], bbSfc[[x]])},
    simplify = FALSE, USE.NAMES = TRUE)
minorRoadsChopped <- sapply(names(bbSfc), function(x) {
    st_intersection(minorRoads[[x]]$osm_lines$geometry, bbSfc[[x]])},
    simplify = FALSE, USE.NAMES = TRUE)
majorRoadsChopped <- sapply(names(bbSfc), function(x) {
    st_intersection(majorRoads[[x]]$osm_lines$geometry, bbSfc[[x]])},
    simplify = FALSE, USE.NAMES = TRUE)
waterChopped <- sapply(names(bbSfc), function(x) {
    a <- st_intersection(water[[x]]$osm_multipolygons$geometry, bbSfc[[x]])
    b <- st_intersection(water[[x]]$osm_polygons$geometry, bbSfc[[x]])
    c <- st_union(a, b)
    return(c)},
    simplify = FALSE, USE.NAMES = TRUE)

## Obtain Red Rooster logo
logo  <- image_read("https://www.redrooster.com.au/favicon.ico")
logo <- image_convert(logo[1],"png")
image_write(logo, "RedRooster.png")

RedRoosterLocations <- sapply(names(bbSfc), function(x) {
    a <- data.frame(st_coordinates(fastfood[[x]]$osm_points[which(fastfood[[x]]$osm_points$brand == "Red Rooster"), ]))
    a$image <- paste0(getwd(),"/RedRooster.png")
    return(a)},
    simplify = FALSE, USE.NAMES = TRUE)

                                  
## Combine it all in a plot for a city, with the Red Rooster data as points

plotCity  <- function(x, logoSize = .02) {
    ggplot() +
        geom_sf(data = majorRoadsChopped[[x]],
                inherit.aes = FALSE,
                color = "grey50",
                size = 0.2) +
        geom_sf(data = boundariesChopped[[x]],
                inherit.aes = FALSE,
                fill = "bisque",
                color = "grey80",
                alpha = .1
                ) +
        geom_sf(data = minorRoadsChopped[[x]],
                inherit.aes = FALSE,
                color = "grey90",
                size = 0.1) +
        geom_sf(data = waterChopped[[x]],
                inherit.aes = FALSE,
                fill = "cadetblue3",
                ##alpha = .1,
                color = NA) +
        geom_image(data = RedRoosterLocations[[x]], aes(x = X, y = Y, image = image), size = logoSize) +
        geom_sf(data = ocean[[x]], fill = "cadetblue3") +
        theme_void() +
        labs(colour = "")
}

## Plot the results and name the CBDs

dev.new(width = 10 * sfext::sf_bbox_asp(bb$Perth), height = 10, unit = "cm")
p <- plotCity("Perth", .015) +
    geom_sf_text(data = boundaries$Perth$osm_multipolygons[which(boundaries$Perth$osm_multipolygons$name == "Perth"), ], label = "Perth", nudge_y = -0.01)
p
png("plots/Perth_Red_Rooster_locations.png", width = 1000 * sfext::sf_bbox_asp(bb$Perth), height = 1000, res = 135)
print(p)
dev.off()
system("mogrify -trim plots/Perth_Red_Rooster_locations.png")
system(paste0("mogrify -resize ", round(1000 * sfext::sf_bbox_asp(bb$Perth), 0), "x1000 plots/Perth_Red_Rooster_locations.png"))

dev.new(width = 10 * sfext::sf_bbox_asp(bb$Sydney), height = 10, unit = "cm")
plotCity("Sydney") + geom_sf_text(data = boundaries$Sydney$osm_multipolygons[which(boundaries$Sydney$osm_multipolygons$name=="Sydney"),], label = "Sydney")

## Cut a line across the city of Sydney, from Windsor in the north-west to Sydney Airport in the south-east, and a curious trend emerges.

Windsor <- st_centroid(boundaries$Sydney$osm_multipolygons[which(boundaries$Sydney$osm_multipolygons$name=="Windsor"),]$geometry)

SydneyAirport <- st_centroid(boundaries$Sydney$osm_multipolygons[which(boundaries$Sydney$osm_multipolygons$name=="Mascot"),]$geometry)

RedRoosterLine <- st_cast(st_union(Windsor,SydneyAirport) ,"LINESTRING")

dev.new(width = 10 * sfext::sf_bbox_asp(bb$Sydney), height = 10, unit = "cm")
p <- plotCity("Sydney") +
    geom_sf(data = RedRoosterLine, color = "#91131B", linetype = "twodash", linewidth = 2) +
    geom_sf_text(data=Windsor, label = "Windsor", nudge_x = -0.015, nudge_y = -0.015) +
    geom_sf_text(data=SydneyAirport, label = "Sydney Airport", nudge_y = -0.005)

p
png("plots/Sydney_Red_Rooster_line.png", width = 800 * sfext::sf_bbox_asp(bb$Sydney), height = 800, res = 135)
print(p)
dev.off()
system("mogrify -trim plots/Sydney_Red_Rooster_line.png")
system(paste0("mogrify -resize ", round(800 * sfext::sf_bbox_asp(bb$Sydney), 0), "x800 plots/Sydney_Red_Rooster_line.png"))