library(sf) library(dplyr) library(ggplot2) library(maptiles) library(tidyterra) library(viridis) library(httr2) ## Obtain the shapefile from WA opendata (SLIP) tempdirSHP <- tempdir() tempfileSHP <- tempfile() req <- request("https://sso.slip.wa.gov.au/as/token.oauth2") |> req_headers("Authorization" = "Basic ZGlyZWN0LWRvd25sb2Fk") |> req_body_form(grant_type = "password", # SLIP username and password stored in # pass - the standard unix password manager username = system2("pass", args = "slip.wa.gov.au | grep Username | sed -e 's/Username: //'", stdout = TRUE), password = system2("pass", args = "slip.wa.gov.au | head -1", stdout = TRUE)) tokenResponse <- req_perform(req) slipUrl <- "https://direct-download.slip.wa.gov.au/datadownload/Education/Current_Active_Schools_Sem_1_2022_Public_DET_020_WA_GDA94_Public_Shapefile.zip" req <- request(slipUrl) |> req_headers( 'Authorization' = paste0('Bearer ',resp_body_json(tokenResponse)$access_token)) responseSlip <- req_perform(req) writeBin(resp_body_raw(responseSlip), con = tempfileSHP) unzip(tempfileSHP, exdir = tempdirSHP) shape.df <- st_read(dsn = tempdirSHP) shape.df$latitude <- as.numeric(shape.df$latitude) shape.df$longitude <- as.numeric(shape.df$longitude) shape.df$totalschoo <- as.numeric(shape.df$totalschoo) shape.df$totalsecon <- as.numeric(shape.df$totalsecon) # filter out remote islands shape.df <- filter(shape.df,longitude > 100) # Obtain the contour box of the shape lat <- c(min(shape.df$latitude), max(shape.df$latitude)) lon <- c(min(shape.df$longitude), max(shape.df$longitude)) # Define Perth metro contours lat.metro <- c(-32.1547, -31.8519) lon.metro <- c(115.7162,116.0733) # Obtain an openstreetmap background image sfg_WA <- st_polygon(list(rbind(c(lon[1],lat[1]), c(lon[2],lat[1]), c(lon[2],lat[2]), c(lon[1],lat[2]), c(lon[1],lat[1]) ))) WATiles <- get_tiles(st_transform(st_sfc(sfg_WA, crs = 4326), 3857), provider = "Stadia.Stamen.TonerLite", crop = TRUE, zoom = 6) ggplot() + geom_spatraster_rgb(data = WATiles) + geom_sf(data = shape.df, aes(color = totalschoo), show.legend = TRUE) + scale_color_viridis(option="viridis", name = "Number\nof pupils") + theme(legend.position = "right", axis.text = element_blank(), axis.ticks = element_blank()) + labs(title = "Western Australian schools", caption = "Map tiles by Stamen Design, under CC BY 3.0.\nMapdata by OpenStreetMap, under ODbL. School data © Government of WA 2022.\n Dep. of Education & Training.", x="", y="") # Zooming in on Perth metro shape.metro.df <- filter(shape.df, abs(latitude) <= max(abs(lat.metro)), abs(latitude) >= min(abs(lat.metro)), abs(longitude) <= max(abs(lon.metro)), abs(longitude) >= min(abs(lon.metro))) sfg_metro <- st_polygon(list(rbind(c(lon.metro[1],lat.metro[1]), c(lon.metro[2],lat.metro[1]), c(lon.metro[2],lat.metro[2]), c(lon.metro[1],lat.metro[2]), c(lon.metro[1],lat.metro[1]) ))) metroTiles <- get_tiles(st_transform(st_sfc(sfg_metro, crs = 4326), 3857), provider = "CartoDB.PositronNoLabels", crop = TRUE, zoom = 11) # Let's look at secondary school only # 95th percentile of school size secon95th <- quantile(shape.metro.df[,"totalsecon">0]$totalsecon,probs = .95) ggplot() + geom_spatraster_rgb(data = metroTiles) + geom_sf(data = filter(shape.metro.df, totalsecon > 0), aes(color = totalsecon, size = totalsecon) ) + geom_sf_text(data = filter(shape.metro.df, totalsecon > secon95th), # name top 5% schools aes(x = longitude, y = latitude, label = schoolname), size = 2) + scale_color_viridis(option="cividis", name = "Number\nof pupils") + guides(size = 'none') + # Do not show legend for size, as color is the main key theme_minimal() + theme(legend.position = "right", panel.grid.major = element_line(colour = "white"), axis.ticks = element_blank()) + labs(title = "Western Australian secondary schools", subtitle = "Perth metro area", caption = "Map tiles © CARTO. Mapdata by OpenStreetMap, under ODbL.\nSchool data © Government of WA, Dep. of Education & Training.", x="", y="") ggsave("WA_schools.png", units = "cm", width = 15, height = 15, dpi = 300, bg = "white", scale = 1.5)