From: Frederik Vanrenterghem <frederik@vanrenterghem.biz>
Date: Tue, 27 Feb 2024 14:05:11 +0000 (+0800)
Subject: Commit analysis code used for blog article.
X-Git-Url: http://git.vanrenterghem.biz/R/openstreetmap-red-rooster-locations-analysis.git/commitdiff_plain/35972914a8824f55bb778f1d3733588b9589f89e

Commit analysis code used for blog article.
---

35972914a8824f55bb778f1d3733588b9589f89e
diff --git a/RedRooster.png b/RedRooster.png
new file mode 100644
index 0000000..8b16927
Binary files /dev/null and b/RedRooster.png differ
diff --git a/openstreetmap_red_rooster_locations_analysis.R b/openstreetmap_red_rooster_locations_analysis.R
new file mode 100644
index 0000000..dc3658d
--- /dev/null
+++ b/openstreetmap_red_rooster_locations_analysis.R
@@ -0,0 +1,184 @@
+## 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"))