openstreetmap_red_rooster_locations_analysis.R

   1 ## The Red Rooster line
   2
   3 ## Initiated by a radio interview described on https://www.abc.net.au/news/2023-11-25/disadvantaged-students-medicine-university/103141050
   4
   5 library(osmdata)
   6 library(sf)
   7 library(sfext) # Used to quickly calculate aspect ratio of bounding box. Install with pak::pkg_install("elipousson/sfext")
   8 library(maptiles)
   9 library(tidyterra)
  10 library(ggplot2)
  11 library(pals)
  12 library(ggimage)
  13 library(magick)
  14
  15 ## Obtain boundaries for significant urban areas in Australia
  16 ## https://www.abs.gov.au/statistics/standards/australian-statistical-geography-standard-asgs-edition-3/jul2021-jun2026/access-and-downloads/digital-boundary-files
  17
  18 shapeAustralia <- read_sf("~/Downloads/SUA_2021_AUST_GDA94.shp")
  19
  20 ## Create bounding boxes around the city shapes
  21 bb <- list()
  22 bb[["Perth"]] <- st_bbox(st_transform(st_as_sf(shapeAustralia[shapeAustralia$SUA_NAME21 == "Perth",]), 4326))
  23 bb[["Sydney"]] <- st_bbox(st_as_sf(shapeAustralia[shapeAustralia$SUA_NAME21 == "Sydney",]))
  24
  25 ## Obtain OpenStreetMap data on the location of fast food restaurants across the cities
  26
  27 fastfood <- lapply(bb, function(x) {
  28     opq(x, timeout = 50) |>
  29     add_osm_feature("amenity", "fast_food") |>
  30     add_osm_feature("access", "!private") |>
  31     osmdata_sf()})
  32
  33 ## Map everything using downloaded background tiles
  34
  35 bbsfPerth <- st_as_sf(shapeAustralia[shapeAustralia$SUA_NAME21 == "Perth",])
  36  TODO - find B/W tile server
  37 mapPerth <- get_tiles(x = bbsfPerth, provider = "OpenStreetMap", zoom = 9)
  38 dev.new(width = 15, height = 15 * 9/16, unit = "cm")
  39 ggplot(fastfood[["Perth"]]$osm_points) +
  40     geom_spatraster_rgb(data = mapPerth) +
  41     geom_sf(aes(colour = brand),
  42             size = 2) +
  43     theme_void() +
  44     scale_color_manual(values=as.vector(glasbey())) +
  45     theme(legend.position="right") +
  46     labs(title = "Fastfood in Perth metro",
  47          colour = "Brand")
  48
  49 ## Obtain way data from OpenStreetMap in order to construct our own background map
  50
  51 majorRoads <- lapply(bb, function(x) {
  52     opq(x, timeout = 120) |>
  53     add_osm_feature(key = "highway",
  54                   value = c("motorway", "primary", "secondary")) |>
  55     osmdata_sf()})
  56
  57 minorRoads <- lapply(bb, function(x) {
  58     opq(x, timeout = 120) |>
  59         add_osm_feature(key = "highway", value = c("tertiary")) |>
  60         osmdata_sf()})
  61
  62 boundaries <- lapply(bb, function(x) {
  63     opq(x, timeout = 120) |>
  64         add_osm_feature(key = "boundary", value = c("administrative")) |>
  65         add_osm_feature(key = "admin_level", value = c(2,8,9,10,11)) |>
  66         osmdata_sf()
  67 })
  68
  69 water <- lapply(bb, function(x) {
  70     opq(x, timeout = 180) |>
  71     add_osm_feature(key = "natural", value = "water") |>
  72     add_osm_feature(key = 'name') |>
  73     osmdata_sf()})
  74
  75 ## Create city outlines, sfc bounding box and separate out ocean area as OSM does not have that polygon
  76
  77 outlines <- sapply(names(bb), function(x) {
  78     st_union(boundaries[[x]]$osm_multipolygons[which(boundaries[[x]]$osm_multipolygons$admin_level == 9), ])},
  79     simplify = FALSE, USE.NAMES = TRUE)
  80 bbSfc <- sapply(names(bb), function(x) {
  81     a <- st_as_sfc(bb[[x]])
  82     a <- st_transform(a, crs = st_crs(outlines[[x]]))
  83     return(a)},
  84     simplify = FALSE, USE.NAMES = TRUE)
  85 ocean <- sapply(names(bbSfc), function(x) {
  86     a <- st_difference(bbSfc[[x]],outlines[[x]])
  87     a <- st_transform(a, crs = st_crs(outlines[[x]]))
  88     return(a)},
  89     simplify = FALSE, USE.NAMES = TRUE)
  90
  91 ## Chop all the osm data to the sfc bounding box, as the queried data runs over its bounding box for many features
  92
  93 boundariesChopped <- sapply(names(bbSfc), function(x) {
  94     st_intersection(boundaries[[x]]$osm_multipolygons[which(boundaries[[x]]$osm_multipolygons$admin_level == 9), ], bbSfc[[x]])},
  95     simplify = FALSE, USE.NAMES = TRUE)
  96 minorRoadsChopped <- sapply(names(bbSfc), function(x) {
  97     st_intersection(minorRoads[[x]]$osm_lines$geometry, bbSfc[[x]])},
  98     simplify = FALSE, USE.NAMES = TRUE)
  99 majorRoadsChopped <- sapply(names(bbSfc), function(x) {
 100     st_intersection(majorRoads[[x]]$osm_lines$geometry, bbSfc[[x]])},
 101     simplify = FALSE, USE.NAMES = TRUE)
 102 waterChopped <- sapply(names(bbSfc), function(x) {
 103     a <- st_intersection(water[[x]]$osm_multipolygons$geometry, bbSfc[[x]])
 104     b <- st_intersection(water[[x]]$osm_polygons$geometry, bbSfc[[x]])
 105     c <- st_union(a, b)
 106     return(c)},
 107     simplify = FALSE, USE.NAMES = TRUE)
 108
 109 ## Obtain Red Rooster logo
 110 logo  <- image_read("https://www.redrooster.com.au/favicon.ico")
 111 logo <- image_convert(logo[1],"png")
 112 image_write(logo, "RedRooster.png")
 113
 114 RedRoosterLocations <- sapply(names(bbSfc), function(x) {
 115     a <- data.frame(st_coordinates(fastfood[[x]]$osm_points[which(fastfood[[x]]$osm_points$brand == "Red Rooster"), ]))
 116     a$image <- paste0(getwd(),"/RedRooster.png")
 117     return(a)},
 118     simplify = FALSE, USE.NAMES = TRUE)
 119
 120
 121 ## Combine it all in a plot for a city, with the Red Rooster data as points
 122
 123 plotCity  <- function(x, logoSize = .02) {
 124     ggplot() +
 125         geom_sf(data = majorRoadsChopped[[x]],
 126                 inherit.aes = FALSE,
 127                 color = "grey50",
 128                 size = 0.2) +
 129         geom_sf(data = boundariesChopped[[x]],
 130                 inherit.aes = FALSE,
 131                 fill = "bisque",
 132                 color = "grey80",
 133                 alpha = .1
 134                 ) +
 135         geom_sf(data = minorRoadsChopped[[x]],
 136                 inherit.aes = FALSE,
 137                 color = "grey90",
 138                 size = 0.1) +
 139         geom_sf(data = waterChopped[[x]],
 140                 inherit.aes = FALSE,
 141                 fill = "cadetblue3",
 142                 ##alpha = .1,
 143                 color = NA) +
 144         geom_image(data = RedRoosterLocations[[x]], aes(x = X, y = Y, image = image), size = logoSize) +
 145         geom_sf(data = ocean[[x]], fill = "cadetblue3") +
 146         theme_void() +
 147         labs(colour = "")
 148 }
 149
 150 ## Plot the results and name the CBDs
 151
 152 dev.new(width = 10 * sfext::sf_bbox_asp(bb$Perth), height = 10, unit = "cm")
 153 p <- plotCity("Perth", .015) +
 154     geom_sf_text(data = boundaries$Perth$osm_multipolygons[which(boundaries$Perth$osm_multipolygons$name == "Perth"), ], label = "Perth", nudge_y = -0.01)
 155 p
 156 png("plots/Perth_Red_Rooster_locations.png", width = 1000 * sfext::sf_bbox_asp(bb$Perth), height = 1000, res = 135)
 157 print(p)
 158 dev.off()
 159 system("mogrify -trim plots/Perth_Red_Rooster_locations.png")
 160 system(paste0("mogrify -resize ", round(1000 * sfext::sf_bbox_asp(bb$Perth), 0), "x1000 plots/Perth_Red_Rooster_locations.png"))
 161
 162 dev.new(width = 10 * sfext::sf_bbox_asp(bb$Sydney), height = 10, unit = "cm")
 163 plotCity("Sydney") + geom_sf_text(data = boundaries$Sydney$osm_multipolygons[which(boundaries$Sydney$osm_multipolygons$name=="Sydney"),], label = "Sydney")
 164
 165 ## 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.
 166
 167 Windsor <- st_centroid(boundaries$Sydney$osm_multipolygons[which(boundaries$Sydney$osm_multipolygons$name=="Windsor"),]$geometry)
 168
 169 SydneyAirport <- st_centroid(boundaries$Sydney$osm_multipolygons[which(boundaries$Sydney$osm_multipolygons$name=="Mascot"),]$geometry)
 170
 171 RedRoosterLine <- st_cast(st_union(Windsor,SydneyAirport) ,"LINESTRING")
 172
 173 dev.new(width = 10 * sfext::sf_bbox_asp(bb$Sydney), height = 10, unit = "cm")
 174 p <- plotCity("Sydney") +
 175     geom_sf(data = RedRoosterLine, color = "#91131B", linetype = "twodash", linewidth = 2) +
 176     geom_sf_text(data=Windsor, label = "Windsor", nudge_x = -0.015, nudge_y = -0.015) +
 177     geom_sf_text(data=SydneyAirport, label = "Sydney Airport", nudge_y = -0.005)
 178
 179 p
 180 png("plots/Sydney_Red_Rooster_line.png", width = 800 * sfext::sf_bbox_asp(bb$Sydney), height = 800, res = 135)
 181 print(p)
 182 dev.off()
 183 system("mogrify -trim plots/Sydney_Red_Rooster_line.png")
 184 system(paste0("mogrify -resize ", round(800 * sfext::sf_bbox_asp(bb$Sydney), 0), "x800 plots/Sydney_Red_Rooster_line.png"))