From 1c447d8b9929265997c51ed9156ca1b187ba84bb Mon Sep 17 00:00:00 2001
From: Thomas Rushby <t.w.rushby@soton.ac.uk>
Date: Wed, 22 Jul 2020 11:44:11 +0100
Subject: [PATCH] Update example to use data from resource drive and tidy/add
annotations.
---
howTo/gisExample.R | 86 +++++++++++++++++++++++++++++-----------------
1 file changed, 54 insertions(+), 32 deletions(-)
diff --git a/howTo/gisExample.R b/howTo/gisExample.R
index fadcbb1..4546c7e 100644
--- a/howTo/gisExample.R
+++ b/howTo/gisExample.R
@@ -1,71 +1,87 @@
-#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
--
GitLab