From 35972914a8824f55bb778f1d3733588b9589f89e Mon Sep 17 00:00:00 2001
From: Frederik Vanrenterghem <frederik@vanrenterghem.biz>
Date: Tue, 27 Feb 2024 22:05:11 +0800
Subject: [PATCH 1/1] Commit analysis code used for blog article.
---
RedRooster.png | Bin 0 -> 2347 bytes
...streetmap_red_rooster_locations_analysis.R | 184 ++++++++++++++++++
2 files changed, 184 insertions(+)
create mode 100644 RedRooster.png
create mode 100644 openstreetmap_red_rooster_locations_analysis.R
diff --git a/RedRooster.png b/RedRooster.png
new file mode 100644
index 0000000000000000000000000000000000000000..8b16927f43557c69f8fe0cacb30fd987d70b449f
GIT binary patch
literal 2347
zcmZ{mdo<JkAICrBzWTaE`Y~l;NLVw`>YL1Ub0-lik!jRyxrHUkM&Gn3lDXv4SRzWb
zax1wcel82S-_~fzrA3+RzWe_3JLmVu_q@*g^*rbOdOlvS^E&7K$NMi=XS-c;YH|Po
z>_XdHyNeX|&ykT5WiRoyWRYlIuy=O?fLJI1+@Ju!y6Ec03;;xmjC1D!0GR~<s-)Z+
zx8tIMRDgq>H6Z*a3hJ5Z03fpuZEbl9KfXXGldzr|tron|+NH6}Rc$4yuiBu&>GxM?
z_xGrk#lz04++!BO9gAN{Zd%!D1i9LR!SpLq(iaud{xk25iyq2+h`j5bBDr6KGK{iE
z5$t<<jua2s{f@l4RNrdBs?@GUvv%3$UKb8llb5O|mN&-3<152C;zpBonvY}$BmsHN
zpy90WlRnZcv0o;Hwxa<jrH%4~4Z3zb3o9`+_rof{dfB>lPnFxN#7qyTD9GqrIzrO)
zsZ|8mbkYHZ;i=qWNj<98-oXB0=elrS1~{msS6lv_OL-ZH+YCe9`kH=_TBTBvt?*!h
zAW+wuIY&f~uW+=$<X^^9yZ@54QC2Nb&K~LPAukoqUNVCq2#w71B__?aHReB3{NRGG
z4uiI#D+~<etkE~WI2BO;1PZrdgWuQP48moDH-CA|*RP%^zszFXDiY6d7Y*opRvYHd
z6^p1$o^!!E;>gETK&tFt>G_Z7r}Po2$t+@AdZv>)dvp=K22d(I-YjCq$pIk#)mo5&
zuS83pw?3lsVtabt&2LwFV_SpN-pBSnS#~mur}mGSWp?wz3S&lhuRF}@F`R`NlCxn|
zB<?eTmSRTmS7JyYC+U|x^X26Vlc<cO6Q-Fa&zED8&TSZc-XiYC_BEEj&PRn;+tOqa
zwoex7?nD!fUm^!4Q>{Ppy@w99o*Q$jV<xkEiYdmvv<UO|elOOI?Tc^t{q?K&Qd`(t
z*E6rJhWZ>CI})3waLOpu3YjO$#&b8MS(m35q_5whJ95h9sy34?3iu_)4-5;+t%Gnc
zSBMjqMtiJgkU{Xa1{M!g7`hXj0Jslakr3xrOmk(WSPKn~P{o=gMcU!#TWdGD@fqTj
zlERH^%sh-fb#eK?3LivI(G~3x#ur5IJ8QKX<HpSO(2uCRicwsuldw-tq^kGtuEQ(#
z%rXKjO+dujNG!^o<+S3rA%coFeyhNj&amT7_#6q<WRIST?i?QWCVgKQ-j(l|!&>n?
ze9^yuam(^u3J$|V9gCO`3kS~G%^|I7{Vq;5yf4L{F%Rh&?LKGOQgpf}G5Lp?>DL^K
z%FEntG=Q2U;u|K|t}9h^=icQRE<;Klf}jNofj?6>b05f$S{&Sv)0ah<DdIzsx=k+r
zK4qV;1yCxa_XkXxjksKNhDdC66S#TO`wv*a;d#Acyee@W(k4+@u)(o2uUYQa>qYKH
zQX9IaXV!);XxEVOTOSuFUXhMDH-dc=1hr$}INTjRQ&XZHar7rMPh=?R6a=d%to_{8
z=XbS}&Ni3A=&Oaky_=+7DV$I)AlJW&=NUcCG_*O>kLMbuuy5mD#cx!|l5eD*JQ?F_
zvu#h&SskU7l}O?89iv#ZW8k2+Kd*8xWo?h;pd~W=XkFVKK0Zr>l)l#}EiQ>2%Yz2R
z#}bCHA#%&As+tnj_zdi!6Z#+C`W6R%7$EwO^9KBV_bgYegFC+y=|1m|ORah5ViBjW
zJLogCO=0#+k5hYNc|o}!Zp+5XwG4VW;nfV^5Stx~@nhH2QJ9}vZSY5^!bt<|<`0cY
zmZgQ=l<$~qD5F_5xWUuYvyqf$q#89YF`}kreb1$YB-l##KDI3e>}_gdtV@6V;^yI#
z;PKok-P*-ikDB%qK~GbXv<|I%RPBZBdbbg8{qy|9nRXX7M&IlO+=z_6B1*ZauWby}
zIo;UBb7P&|`#6|cb`~<tPzgbcp)fa}8AkKI_?PrF*VM>CdJsqXA9F-+|Ke$$TF3zQ
z_Tn<9e5aDn3vmv+?HWnqzc3j?UrdDX*x=E`xfCy@rCFWuDH)EZMizZ%dsXTCZy@ks
zmC6Y037)`dwwP#Ek>k)!q)i?K2%gmsR%7Yn8I5jkA9|MJQ?u{7+)rsSUvD=O?{s>5
zscGkYbi(?j_qT9+sh%$gAxPxgDNeFNw)XER)xjgN4XJ$^S*6N7vnTCoGuvBzQ5BHS
zV~LBa9<MLUek5f)ZCMbt2)w(CQ~pGAkR@}$y#10{rSUu_<459sw@<fh(JR02tL>ed
z?AI+l6DeDf9rAT4%429eX>l1Y*i3m+C##H+HK7Oopgy!;SyA$4cN9L1Pkw(Gg`W=l
zqkfY-`r0)5`y9KT8U!V8nwM^iRav#lKEpak&6qb47jTTq)WC@iaCCqJM9Eo!im7=6
zP)Z&jN_of3*NGFrGX2@B9pMal4cq1DA070>h-%X5(+3zQXOL})ap27_>6$u=*=psE
zQV^ZRVOS8<Aqs!X%2=4mzCA7Y&eO_s9yR@9&rH`kc9qYWSk+`ScTGnNXL+AmxjL&9
z;FbC{l~t8E>P^eJ%6c^yf+4+Wzru*{kJk??8YlWbNyiJc#;e1BbUdQbF-l{#@^&!$
zh2k=jggAJt8U)ede9-`tBhB_N4q3hBDml-`Bw8<#DPENaK5y)6jan$FA}DIb#d((^
zL`LtVgvvrOFFgv~`LW|}t{_=wH3^gC2GQwydXQ4#Uu-g!j>dk0%-5k$>@heVdA6x5
zC#y}nXBs)wU=Y^irXVgt6=_9>KweI!YtdkwkGYj+*6hyHP@2IhrxW+|Rz|Mxul5Gp
zX}V^a<+OCEc_+DrNkThz=tV*?`J>RldHxxp9j(|k;r7c^2fp%Av&=GOXegetH!q+I
zYT&+{a^+1*QL9bN$K0F{4*!L)fgcZPFehfx6AO3btO>5`F6IaCFq*uODg=<~y0v4P
zE}^2ogHD8vcSN9n1P*y23?~u*4uc~MVQ@qEkyD2eNW@Vj%)|f&L&9K*tR|KJ0|>;x
s;Gmd)2Y!Ahmx=(r|4oPpCg8#&{0XFgXNtFdqC^aUwsE$uK%Kw&H`qvIO#lD@
literal 0
HcmV?d00001
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"))
--
2.39.5