+## 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"))
