Update example to use data from resource drive and tidy/add annotations.

howTo/gisExample.R

-#Install/load pacman
-if(!require(pacman)){install.packages("pacman");require(pacman)}
-#Install/load tons of packages
-p_load(
-  leaflet,
-  sf,
-  ggplot,
-  data.table,
-  ggspatial,
-  dplyr
+#### An example of plotting Census Output areas in R with ggplot and leaflet
+#### created by Tom Rushby (t.w.rushby@soton.ac.uk, @tom_rushby)
+
+# Libraries ----
+library(woRkflow) # remember to build it first :-)
+woRkflow::setup() # load env.R set up the default paths etc
+
+# list required packages
+reqLibs <- c("leaflet",
+             "sf",
+             "ggplot2",
+             "data.table",
+             "ggspatial",
+             "dplyr",
+             "drake", # what's done stays done
+             "here", # here
+             "lubridate", # dates and times
+             "bookdown" # for making reports (should also install knittr etc)
 )
+# install/load them
+woRkflow::loadLibraries(reqLibs)
 
+# Set data file paths ----
 
+# Set path for geography data
+dPathGeo <- paste0(dPath,"geography/census/")
+# Set path for SAVE area-based estimated electricity demand profile data
+dPathSAVE <- sub("/data", "/SAVE/data", dPath,fixed = TRUE)
 
+# Load data ----
 
-# Census 2011 ----
+# Census 2011
 # https://www.cultureofinsight.com/post/multivariate-dot-density-maps-in-r-with-sf-ggplot2
 
-# Loading 2011 Census output area shape file
+# Load 2011 Census output area shape file
 # Downloaded from https://geoportal.statistics.gov.uk/
 
-
-dPathGeo <- "/Users/twr1m15/documents/data/geographical/census/"
-# ukOA2011 <- st_read(paste0(dPathGeo,"England_oac_2011_shp/england_oac_2011.shp"),stringsAsFactors = FALSE,quiet = TRUE)
-
 # use 'clipped' boundaries - these exclude areas of water
 ukOA2011b <- st_read(paste0(dPathGeo,"Output_Area_December_2011_Generalised_Clipped_Boundaries_in_England_and_Wales/Output_Area_December_2011_Generalised_Clipped_Boundaries_in_England_and_Wales.shp"),stringsAsFactors = FALSE,quiet = TRUE)
 
-# Load Census area lookup tables
+# Load Census area lookup tables for matching to regions
 ukLookupOA <- read.csv(paste0(dPathGeo,"Lookup/PCD_OA_LSOA_MSOA_LAD_MAY19_UK_LU.csv"), header=TRUE, sep=",", na.strings= c("NA", " ", ""), stringsAsFactors = FALSE)
 ukLookupReg <- read.csv(paste0(dPathGeo,"Lookup/Output_Area_to_Region_December_2018_Lookup_in_England_and_Wales.csv"), header=TRUE, sep=",", na.strings= c("NA", " ", ""), stringsAsFactors = FALSE)
 
 head(ukLookupOA)
 
+# remove a bunch of columns we don't need
 ukLookupOA <- ukLookupOA %>%
   select(-c("pcd7","pcd8","pcds","dointr","doterm","usertype","ladnmw"))
 
+# filter out NA regions and remove unwanted columns
 ukLookupReg <- ukLookupReg %>%
   dplyr::filter(!is.na(RGN18NM)) %>%
   select(-c("RGN18NMW", "FID"))
 
-
+# remove duplicates
 ukLookupOA <- unique(ukLookupOA)
 
+# 232,038 unique OAs
 uniqueN(ukLookupOA$oa11cd)
 
 head(ukLookupReg)
 
-# join OA lookup data and filter by region and local authority
+# join OA lookup data with region lookup
 ukLookup <- left_join(ukLookupOA,ukLookupReg, by = c("oa11cd" = "OA11CD"))
-
+# and filter by region (South East) and local authority (Southampton)
 ukLookup <- ukLookup %>%
   dplyr::filter(RGN18NM == "South East") %>%
   dplyr::filter(ladnm == "Southampton")
 
 head(ukLookup)
+uniqueN(ukLookup$oa11cd) # now 766 OAs.
 
 # merge OA shape and lookup data
 sf_data <- left_join(ukLookup,ukOA2011b, by = "oa11cd") %>%
   st_as_sf() # reset as sf class
 
 # remove data tables no longer req'd
-#rm(list = ls(pattern = "^ukLookup"))
-#rm(ukOA2011)
+rm(list = ls(pattern = "^ukLookup"))
+rm(ukOA2011b)
 
-# Plot the map ----
+# Plot the map (ggplot) ----
 
-# Using ggplot2 ----
+# Using ggplot2
 # see https://ggplot2.tidyverse.org/reference/ggsf.html
 # some features like scale bar and north point require ggspatial package
 
@@ -77,7 +93,7 @@ ggplot(sf_data) +
 
 # Make interactive through plotly? https://plot.ly/ggplot2/maps-sf/
 
-# Using leaflet ----
+# Plot the map (leaflet) ----
 # Useful help for leaflet
 # https://rstudio.github.io/leaflet/choropleths.html
 library(htmltools)
@@ -89,14 +105,14 @@ st_crs(sf_data) # current coord system EPSG: 27700
 # Useful lookup spatial reference for CRS
 # https://spatialreference.org/ref/epsg/27700/
 
-st_crs(sf_data) # now CRS shows EPSG: 4326 (what leaflet wants)
-
-# transform the
+# transform the coord system
 sf_data <- st_transform(sf_data, "+proj=longlat +datum=WGS84")
 
-# Merge SAVE oa level data ----
-# Load file - created by gisDataProcess.Rmd
-saveOAdata <- read.csv(paste0(dPath,"output/GIS/OAstats/gisExportRv2.csv"),stringsAsFactors = FALSE)
+st_crs(sf_data) # now CRS shows EPSG: 4326 (what leaflet wants)
+
+# Merge SAVE oa level demand estimates data
+# Load file - created by gisDataProcess.Rmd in SAVE project
+saveOAdata <- read.csv(paste0(dPathSAVE,"output/GIS/OAstats/gisExportRv2.csv"),stringsAsFactors = FALSE)
 
 # Join with output area data
 saveOAprofiles <- left_join(saveOAdata,sf_data, by = c("oaCode" = "oa11cd")) %>%
@@ -115,6 +131,7 @@ sf_data2$popup_text <-
         '<br/>', 'Local authority: ', '<b>', sf_data2$ladnm, '</b>', ' ') %>%
   lapply(htmltools::HTML)
 
+# create map using leaflet
 leaflet(sf_data2) %>%
   addTiles() %>%  # Add default OpenStreetMap map tiles
   addPolygons(color = "blue", fillColor = "blue", fillOpacity = 0.2, weight = 1.5, popup = ~(oa11cd), # popups clicked
@@ -131,6 +148,10 @@ leaflet(sf_data2) %>%
 
 dev.off()
 
+# Works pretty well
+
+# Update with more info on popups
+
 # create popup first (requires htmltools)
 saveOAprofiles$popup_text <-
   paste("Output area code: ","<b>", saveOAprofiles$oa11cd, "</b>",
@@ -139,10 +160,11 @@ saveOAprofiles$popup_text <-
         '<br/>', 'Peak hours demand (Household): ', '<b>', round(saveOAprofiles$mean4to8HH,3), '</b>',' kW') %>%
   lapply(htmltools::HTML)
 
+# create shading scale
 # pal = colorQuantile("Reds", domain = saveOAprofiles$mean4to8, 5)
 pal = colorBin("Reds", domain = saveOAprofiles$mean4to8, 5, pretty = TRUE)
 
-
+# redraw map
 leaflet(saveOAprofiles) %>%
   addTiles() %>%  # Add default OpenStreetMap map tiles
   addPolygons(color = "blue", fillColor = ~pal(mean4to8), fillOpacity = 0.9, weight = 1, popup = ~(oa11cd), # popups clicked