From: Frederik Vanrenterghem 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"))