]> git.vanrenterghem.biz Git - R/project-wa-schools.git/blob - WA_active_schools.R
18816aab1f3537dd90e100aa5d23974f19456503
[R/project-wa-schools.git] / WA_active_schools.R
1 library(sf)
2 library(dplyr)
3 library(ggplot2)
4 library(OpenStreetMap)
5 library(viridis)
6 library(httr2)
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)))
19 tokenResponse <- req_perform(req)
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"
23 req <- request(slipUrl) |>
24     req_headers( 'Authorization' = paste0('Bearer ',resp_body_json(tokenResponse)$access_token))
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)
35 # filter out remote islands
36 shape.df <- filter(shape.df,longitude > 100)
37 # 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))
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)
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")
57 mapLatLon <- openproj(map)
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="")
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="")