WA_active_schools.R

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