library(sf)
library(dplyr)
library(ggplot2)
library(OpenStreetMap)
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(list("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))

lat.metro <- c(-32.1547, -31.8519)
lon.metro <- c(115.7162,116.0733)
upperLeft.metro <- c(-31.8519,115.7162) 
lowerRight.metro <- c(-32.1547,116.0733)

# Obtain an openstreetmap background image
map <- openmap(upperLeft = c(lat[2],lon[1]), 
               lowerRight = c(lat[1],lon[2]), 
               minNumTiles = 6,
               type="https://stamen-tiles.a.ssl.fastly.net/toner-lite/{z}/{x}/{y}.png")
map.metro <- openmap(upperLeft = upperLeft.metro,
                     lowerRight = lowerRight.metro,
                     minNumTiles = 6,
                     type="https://stamen-tiles.a.ssl.fastly.net/toner-lite/{z}/{x}/{y}.png")

mapLatLon <- openproj(map)

OpenStreetMap::autoplot.OpenStreetMap(mapLatLon) +
  geom_point(data = shape.df, 
             aes(x = longitude, 
                 y = latitude,
                 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 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)))
mapLatLon.metro <- openproj(map.metro)
# Let's look at secondary school only
# 95th percentile of school size
secon95th <- quantile(shape.metro.df[,"totalsecon">0]$totalsecon,probs = .95)
OpenStreetMap::autoplot.OpenStreetMap(mapLatLon.metro) +
  geom_point(data = filter(shape.metro.df, totalsecon > 0), # geom_sf is not yet in ggplot2 on CRAN, so sticking to this for now
             aes(x = longitude, 
                 y = latitude,
                 color = totalsecon,
                 size = totalsecon)
             ) +
  geom_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="viridis",
                      name = "Number\nof pupils") +
  guides(size = 'none') + # Do not show legend for size, as color is the main key
  theme(legend.position = "right",
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  labs(title = "Western Australian secondary schools",
       subtitle = "Perth metro area",
       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="")