WA_active_schools.R

   1 # install libgdal-dev on debian first
   2 # install libudunits2-dev on debian first
   3 #install.packages("sf")
   4 library(sf)
   5 library(tidyverse)
   6 library(broom)
   7 #install.packages("OpenStreetMap")
   8 library(OpenStreetMap)
   9 library(viridis)
  10
  11 # Read the shapefile obtained from WA opendata (SLIP)
  12 #shape <- readOGR(dsn = "data/", layer = "CurrentActiveSchoolsSemester12017_PublicDET_014")
  13 shape.df <- st_read(dsn = "data/", layer = "CurrentActiveSchoolsSemester12017_PublicDET_014", stringsAsFactors = FALSE)
  14 shape.df$latitude <- as.numeric(shape.df$latitude)
  15 shape.df$longitude <- as.numeric(shape.df$longitude)
  16 shape.df$totalschoo <- as.numeric(shape.df$totalschoo)
  17 shape.df$totalsecon <- as.numeric(shape.df$totalsecon)
  18
  19 # filter out remote islands
  20 shape.df <- filter(shape.df,longitude > 100)
  21 #shape.df <- filter(shape.df, totalschoo < 500)
  22 # Obtain the contour box of the shape
  23
  24 lat <- c(min(shape.df$latitude), max(shape.df$latitude))
  25 lon <- c(min(shape.df$longitude), max(shape.df$longitude))
  26
  27 lat.metro <- c(-32.1547, -31.8519)
  28 lon.metro <- c(115.7162,116.0733)
  29 upperLeft.metro <- c(-31.8519,115.7162)
  30 lowerRight.metro <- c(-32.1547,116.0733)
  31
  32 # Obtain an openstreetmap background image
  33 map <- openmap(upperLeft = c(lat[2],lon[1]),
  34                lowerRight = c(lat[1],lon[2]),
  35                minNumTiles = 6,
  36                type="https://stamen-tiles.a.ssl.fastly.net/toner-lite/{z}/{x}/{y}.png")
  37 map.metro <- openmap(upperLeft = upperLeft.metro,
  38                      lowerRight = lowerRight.metro,
  39                      minNumTiles = 6,
  40                      type="https://stamen-tiles.a.ssl.fastly.net/toner-lite/{z}/{x}/{y}.png")
  41
  42 mapLatLon <- openproj(map)
  43
  44 autoplot(mapLatLon) +
  45   geom_point(data = shape.df,
  46              aes(x = longitude,
  47                  y = latitude,
  48                  color = totalschoo),
  49              show.legend = TRUE) +
  50   scale_color_viridis(option="viridis",
  51                       name = "Number\nof pupils") +
  52   theme(legend.position = "right",
  53         axis.text = element_blank(),
  54         axis.ticks = element_blank()) +
  55   labs(title = "Western Australian schools",
  56        caption = "Map tiles by Stamen Design, under CC BY 3.0.\nMapdata by OpenStreetMap, under ODbL.
  57        School data © Government of WA 2017.\n Dep. of Education & Training.",
  58        x="",
  59        y="")
  60
  61
  62 # Zooming in on metro
  63 shape.metro.df <- filter(shape.df, abs(latitude) <= max(abs(lat.metro)),
  64                          abs(latitude) >= min(abs(lat.metro)),
  65                          abs(longitude) <= max(abs(lon.metro)),
  66                          abs(longitude) >= min(abs(lon.metro)))
  67 mapLatLon.metro <- openproj(map.metro)
  68 # Let's look at secondary school only
  69 # 95th percentile of school size
  70 secon95th <- quantile(shape.metro.df[,"totalsecon">0]$totalsecon,probs = .95)
  71 autoplot(mapLatLon.metro) +
  72   geom_point(data = filter(shape.metro.df, totalsecon > 0), # geom_sf is not yet in ggplot2 on CRAN, so sticking to this for now
  73              aes(x = longitude,
  74                  y = latitude,
  75                  color = totalsecon,
  76                  size = totalsecon)
  77              ) +
  78   geom_text(data = filter(shape.metro.df, totalsecon > secon95th), # name top 5% schools
  79              aes(x = longitude,
  80                  y = latitude,
  81                  label = schoolname),
  82             size = 2) +
  83   scale_color_viridis(option="viridis",
  84                       name = "Number\nof pupils") +
  85   guides(size = FALSE) + # Do not show legend for size, as color is the main key
  86   theme(legend.position = "right",
  87         axis.text = element_blank(),
  88         axis.ticks = element_blank()) +
  89   labs(title = "Western Australian secondary schools",
  90        subtitle = "Perth metro area",
  91        caption = "Map tiles by Stamen Design, under CC BY 3.0.\nMapdata by OpenStreetMap, under ODbL.
  92        School data © Government of WA 2017.\n Dep. of Education & Training.",
  93        x="",
  94        y="")