Some mapping tests using sp, ggplot, coordinate transformation, subsetting polygons...

mapping_tests.R

library(ggmap)               # Changes SpatialPolygonsDataFrame into a data frame
library(maptools)            # Reads the shapefile
library(sp)                  # Assigns a CRS to the shapefile
library(tidyverse)

# Shapefile from http://geoportal.statistics.gov.uk/datasets/afcc88affe5f450e9c03970b237a7999_3

shpfile <- paste0("Datasets/Shapefile/Wards_December_2016_Super_Generalised_",
                  "Clipped_Boundaries_in_Great_Britain.shp")
wards <- readShapeSpatial(shpfile)
proj4string(wards) <- CRS("+init=epsg:27700") # From http://spatialreference.org

# readOGR from the rgdal library imports a coordinate reference system
# if there is one embedded in the shapefile. You can check by proj4string(wards)

# fortify() changes SpatialPolygonsDataFrame into a data frame
# wards <- fortify(wards)

ggplot(wards, aes(long, lat)) +
  geom_polygon(aes(group = group), 
               color = "black", fill = NA) +              # This creates an outline map
  coord_quickmap() +
  theme_bw() +                                            # Clean theme with border
  labs(title = "Midlands Repairs",                        # labs() specifies labels
       subtitle = "Average Yearly Repairs per Ward")

# The above code creates an outline map. Now we can create a chloropleth
# where the wards are trimmed to our data points.

# Ward data from http://neighbourhood.statistics.gov.uk/dissemination/instanceSelection.do?
# JSAllowed=true&Function=&$ph=60_61&CurrentPageId=61&step=2&datasetFamilyId=266&
# instanceSelection=121427&Next.x=23&Next.y=9&nsjs=true&nsck=false&nssvg=false&nswid=1920

ward_data <- read.csv("Datasets/WardData/D120301_1153_2003Admin_WARD.CSV",
                      stringsAsFactors = FALSE)

# One way to filter is to use subset directly on the polygon.
wards_filter <- wards %>% subset(lad16nm %in% c("Derby", "Nottingham", "Erewash", "Broxtowe",
                                                "Gedling", "Rushcliffe", "South Derbyshire",
                                                "Amber Valley", "Ashfield"))

ggplot(wards_filter, aes(long, lat)) +
  geom_polygon(aes(group = group), 
               color = "black", fill = NA) +              
  coord_quickmap() +
  theme_bw() +                                            
  labs(title = "Midlands Repairs",                        
       subtitle = "Average Yearly Repairs per Ward")

# Another way to filter is by a selection of points
notm_points <- data.frame(town = c("Keyworth", "Fiskerton", "Alfreton", "Sherwood"),
                          long = c(-1.0903, -0.4281, -1.3823, -1.1540),
                          lat = c(52.8709, 53.2345, 53.0975, 52.9832))

# Changes notm_points into a SpatialPointsDataFrame
coordinates(notm_points) <- ~ long + lat 

# The coordinate reference system must be the same for both
proj4string(notm_points) <- CRS("+init=epsg:4326")
wards <- spTransform(wards, CRS("+init=epsg:4326"))  # Requires gdal libs to be installed on Linux

# over() returns, for each row of points, the corresponding polygon.
# Might need to cbind back to the points if you want the original data.
ward_intersects <- over(notm_points, wards)
wards_intersected <- wards %>% subset(wd16nm %in% ward_intersects$wd16nm)

ggplot(wards_intersected, aes(long, lat)) +
  geom_polygon(aes(group = group), 
               color = "black", fill = NA) +              
  coord_quickmap() +
  theme_bw() +                                            
  labs(title = "Midlands Repairs",                        
       subtitle = "Average Yearly Repairs per Ward")        # This doesn't look right...

View(tbl_df(wards_intersected))     # Damn! Some wards have the same name. Have to use the codes.

wards_intersected_codes <- wards %>% subset(wd16cd %in% ward_intersects$wd16cd)

ggplot(wards_intersected_codes, aes(long, lat)) +
  geom_polygon(aes(group = group), 
               color = "black", fill = NA) +              
  coord_quickmap() +
  theme_bw() +                                            
  labs(title = "Midlands Repairs",                        
       subtitle = "Average Yearly Repairs per Ward") 

# Area metadata can be left joined onto the spatial data
# and for this we need them both to be a data frame.
# Note that the ID of the spatial data must match the id of the frame from which
# we take the values. So there might be some data wrangling involved.
wards_filter_df <- fortify(wards_filter)
ward_ids <- data.frame(ward = wards_filter$wd16cd, id = wards_filter$objectid)

ward_ids$id <- as.character(ward_ids$id)
ward_ids$ward <- as.character(ward_ids$ward)

# Oh no! The ONS must have changed their ward references!
# Now we have to find some other data with the right references.
merge(ward_data, ward_ids, by.x = "WARD_CODE", by.y = "ward")