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
  18 # filter out remote islands
  19 shape.df <- filter(shape.df,longitude > 100)
  20 #shape.df <- filter(shape.df, totalschoo < 500)
  21 # Obtain the contour box of the shape
  22
  23 lat <- c(min(shape.df$latitude), max(shape.df$latitude))
  24 lon <- c(min(shape.df$longitude), max(shape.df$longitude))
  25
  26 lat.metro <- c(-32.1547, -31.8519)
  27 lon.metro <- c(115.7162,116.0733)
  28 upperLeft.metro <- c(-31.8519,115.7162)
  29 lowerRight.metro <- c(-32.1547,116.0733)
  30
  31 # Obtain an openstreetmap background image
  32 map <- openmap(upperLeft = c(lat[2],lon[1]),
  33                lowerRight = c(lat[1],lon[2]),
  34                minNumTiles = 6,
  35                type="https://stamen-tiles.a.ssl.fastly.net/toner-lite/{z}/{x}/{y}.png")
  36 map.metro <- openmap(upperLeft = upperLeft.metro,
  37                      lowerRight = lowerRight.metro,
  38                      minNumTiles = 6,
  39                      type="https://stamen-tiles.a.ssl.fastly.net/toner-lite/{z}/{x}/{y}.png")
  40
  41 mapLatLon <- openproj(map)
  42
  43 autoplot(mapLatLon) +
  44   geom_point(data = shape.df,
  45              aes(x = longitude,
  46                  y = latitude,
  47                  color = totalschoo),
  48              show.legend = TRUE) +
  49   scale_color_viridis(option="viridis",
  50                       name = "Number\nof pupils") +
  51   theme(legend.position = "right",
  52         axis.text = element_blank(),
  53         axis.ticks = element_blank()) +
  54   labs(title = "Western Australian schools",
  55        caption = "Map tiles by Stamen Design, under CC BY 3.0.\nMapdata by OpenStreetMap, under ODbL.
  56        School data © Government of WA 2017.\n Dep. of Education & Training.",
  57        x="",
  58        y="")
  59
  60
  61 # Zooming in on metro
  62 shape.metro.df <- filter(shape.df, abs(latitude) <= max(abs(lat.metro)),
  63                          abs(latitude) >= min(abs(lat.metro)),
  64                          abs(longitude) <= max(abs(lon.metro)),
  65                          abs(longitude) >= min(abs(lon.metro)))
  66 mapLatLon.metro <- openproj(map.metro)
  67 autoplot(mapLatLon.metro) +
  68   geom_point(data = shape.metro.df, # geom_sf is not yet in ggplot2 on CRAN, so sticking to this for now
  69              aes(x = longitude,
  70                  y = latitude,
  71                  color = totalschoo),
  72              show.legend = TRUE) +
  73   scale_color_viridis(option="viridis",
  74                       name = "Number\nof pupils") +
  75   theme(legend.position = "right",
  76         axis.text = element_blank(),
  77         axis.ticks = element_blank()) +
  78   labs(title = "Western Australian schools",
  79        subtitle = "Perth metro area only",
  80        caption = "Map tiles by Stamen Design, under CC BY 3.0.\nMapdata by OpenStreetMap, under ODbL.
  81        School data © Government of WA 2017.\n Dep. of Education & Training.",
  82        x="",
  83        y="")