WA_active_schools.R

   1 library(sf)
   2 library(dplyr)
   3 library(ggplot2)
   4 library(maptiles)
   5 library(tidyterra)
   6 library(viridis)
   7 library(httr2)
   8
   9 ## Obtain the shapefile from WA opendata (SLIP)
  10 tempdirSHP <- tempdir()
  11 tempfileSHP <- tempfile()
  12 req <- request("https://sso.slip.wa.gov.au/as/token.oauth2") |>
  13     req_headers("Authorization" = "Basic ZGlyZWN0LWRvd25sb2Fk") |>
  14     req_body_form(grant_type = "password",
  15                   # SLIP username and password stored in
  16                   # pass - the standard unix password manager
  17                   username = system2("pass", args = "slip.wa.gov.au | grep Username | sed -e 's/Username: //'", stdout = TRUE),
  18                   password = system2("pass", args = "slip.wa.gov.au | head -1", stdout = TRUE))
  19
  20 tokenResponse <- req_perform(req)
  21
  22 slipUrl <-  "https://direct-download.slip.wa.gov.au/datadownload/Education/Current_Active_Schools_Sem_1_2022_Public_DET_020_WA_GDA94_Public_Shapefile.zip"
  23
  24 req <- request(slipUrl) |>
  25     req_headers( 'Authorization' = paste0('Bearer ',resp_body_json(tokenResponse)$access_token))
  26
  27 responseSlip <- req_perform(req)
  28 writeBin(resp_body_raw(responseSlip), con = tempfileSHP)
  29 unzip(tempfileSHP, exdir = tempdirSHP)
  30 shape.df <- st_read(dsn = tempdirSHP)
  31 shape.df$latitude <- as.numeric(shape.df$latitude)
  32 shape.df$longitude <- as.numeric(shape.df$longitude)
  33 shape.df$totalschoo <- as.numeric(shape.df$totalschoo)
  34 shape.df$totalsecon <- as.numeric(shape.df$totalsecon)
  35
  36 # filter out remote islands
  37 shape.df <- filter(shape.df,longitude > 100)
  38 # Obtain the contour box of the shape
  39 lat <- c(min(shape.df$latitude), max(shape.df$latitude))
  40 lon <- c(min(shape.df$longitude), max(shape.df$longitude))
  41 # Define Perth metro contours
  42 lat.metro <- c(-32.1547, -31.8519)
  43 lon.metro <- c(115.7162,116.0733)
  44 # Obtain an openstreetmap background image
  45 sfg_WA <- st_polygon(list(rbind(c(lon[1],lat[1]),
  46                                    c(lon[2],lat[1]),
  47                                    c(lon[2],lat[2]),
  48                                    c(lon[1],lat[2]),
  49                                    c(lon[1],lat[1])
  50                                    )))
  51 WATiles <- get_tiles(st_transform(st_sfc(sfg_WA, crs = 4326), 3857), provider = "Stadia.Stamen.TonerLite", crop = TRUE, zoom = 6)
  52
  53 ggplot() +
  54   geom_spatraster_rgb(data = WATiles) +
  55   geom_sf(data = shape.df,
  56              aes(color = totalschoo),
  57              show.legend = TRUE) +
  58   scale_color_viridis(option="viridis",
  59                       name = "Number\nof pupils") +
  60   theme(legend.position = "right",
  61         axis.text = element_blank(),
  62         axis.ticks = element_blank()) +
  63   labs(title = "Western Australian schools",
  64        caption = "Map tiles by Stamen Design, under CC BY 3.0.\nMapdata by OpenStreetMap, under ODbL.
  65        School data © Government of WA 2022.\n Dep. of Education & Training.",
  66        x="",
  67        y="")
  68
  69
  70 # Zooming in on Perth metro
  71 shape.metro.df <- filter(shape.df, abs(latitude) <= max(abs(lat.metro)),
  72                          abs(latitude) >= min(abs(lat.metro)),
  73                          abs(longitude) <= max(abs(lon.metro)),
  74                          abs(longitude) >= min(abs(lon.metro)))
  75 sfg_metro <- st_polygon(list(rbind(c(lon.metro[1],lat.metro[1]),
  76                                    c(lon.metro[2],lat.metro[1]),
  77                                    c(lon.metro[2],lat.metro[2]),
  78                                    c(lon.metro[1],lat.metro[2]),
  79                                    c(lon.metro[1],lat.metro[1])
  80                                    )))
  81 metroTiles <- get_tiles(st_transform(st_sfc(sfg_metro, crs = 4326), 3857), provider = "CartoDB.PositronNoLabels", crop = TRUE, zoom = 11)
  82 # Let's look at secondary school only
  83 # 95th percentile of school size
  84 secon95th <- quantile(shape.metro.df[,"totalsecon">0]$totalsecon,probs = .95)
  85 ggplot() +
  86     geom_spatraster_rgb(data = metroTiles) +
  87     geom_sf(data = filter(shape.metro.df, totalsecon > 0),
  88             aes(color = totalsecon,
  89                 size = totalsecon)
  90             ) +
  91     geom_sf_text(data = filter(shape.metro.df, totalsecon > secon95th), # name top 5% schools
  92                  aes(x = longitude,
  93                      y = latitude,
  94                      label = schoolname),
  95                  size = 2) +
  96     scale_color_viridis(option="cividis",
  97                         name = "Number\nof pupils") +
  98     guides(size = 'none') + # Do not show legend for size, as color is the main key
  99     theme_minimal() +
 100     theme(legend.position = "right",
 101           panel.grid.major = element_line(colour = "white"),
 102           axis.ticks = element_blank()) +
 103     labs(title = "Western Australian secondary schools",
 104          subtitle = "Perth metro area",
 105          caption = "Map tiles © CARTO. Mapdata by OpenStreetMap, under ODbL.\nSchool data © Government of WA, Dep. of Education & Training.",
 106        x="",
 107        y="")
 108
 109 ggsave("WA_schools.png", units = "cm", width = 15, height = 15, dpi = 300, bg = "white", scale = 1.5)