From 35972914a8824f55bb778f1d3733588b9589f89e Mon Sep 17 00:00:00 2001 From: Frederik Vanrenterghem 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{mdoYL1Ub0-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~H6Z*a3hJ5Z03fpuZEbl9KfXXGldzr|tron|+NH6}Rc$4yuiBu&>GxM? z_xGrk#lz04++!BO9gAN{Zd%!D1i9LR!SpLq(iaud{xk25iyq2+h`j5bBDr6KGK{iE z5$t<lPnFxN#7qyTD9GqrIzrO) zsZ|8mbkYHZ;i=qWNj<98-oXB0=elrS1~{msS6lv_OL-ZH+YCe9`kH=_TBTBvt?*!h zAW+wuIY&f~uW+=$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{Ppy@w99o*Q$jVCI})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#ur5IJ8QKXn? 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)0ahqYKH 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`~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~LBa9ed 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^ZRVOS8P^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.2