# install libgdal-dev on debian first #install.packages("rgdal") library(rgdal) library(tidyverse) #install.packages("OpenStreetMap") library(OpenStreetMap) library(viridis) # Read the shapefile obtained from WA opendata (SLIP) shape <- readOGR(dsn = "data/", layer = "CurrentActiveSchoolsSemester12017_PublicDET_014") shape.df <- tidy(shape) shape.df$latitude <- as.numeric(as.character(shape.df$latitude)) shape.df$longitude <- as.numeric(as.character(shape.df$longitude)) shape.df$totalschoo <- as.numeric(as.character(shape.df$totalschoo)) # filter out remote islands shape.df <- filter(shape.df,longitude > 100) #shape.df <- filter(shape.df, totalschoo < 500) # 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) autoplot(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 2017.\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) autoplot(mapLatLon.metro) + geom_point(data = shape.metro.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", subtitle = "Perth metro area only", caption = "Map tiles by Stamen Design, under CC BY 3.0.\nMapdata by OpenStreetMap, under ODbL. School data © Government of WA 2017.\n Dep. of Education & Training.", x="", y="")