Skip to content
Snippets Groups Projects
Commit 4356b713 authored by Tom Rushby's avatar Tom Rushby
Browse files

Import R script for gisExample. WIP

parent 7b6063ed
No related branches found
No related tags found
No related merge requests found
#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
)
# Census 2011 ----
# https://www.cultureofinsight.com/post/multivariate-dot-density-maps-in-r-with-sf-ggplot2
# Loading 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
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)
ukLookupOA <- ukLookupOA %>%
select(-c("pcd7","pcd8","pcds","dointr","doterm","usertype","ladnmw"))
ukLookupReg <- ukLookupReg %>%
dplyr::filter(!is.na(RGN18NM)) %>%
select(-c("RGN18NMW", "FID"))
ukLookupOA <- unique(ukLookupOA)
uniqueN(ukLookupOA$oa11cd)
head(ukLookupReg)
# join OA lookup data and filter by region and local authority
ukLookup <- left_join(ukLookupOA,ukLookupReg, by = c("oa11cd" = "OA11CD"))
ukLookup <- ukLookup %>%
dplyr::filter(RGN18NM == "South East") %>%
dplyr::filter(ladnm == "Southampton")
head(ukLookup)
# 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)
# Plot the map ----
# Using ggplot2 ----
# see https://ggplot2.tidyverse.org/reference/ggsf.html
# some features like scale bar and north point require ggspatial package
ggplot(sf_data) +
geom_sf() +
annotation_scale(location = "bl", width_hint = 0.5) +
annotation_north_arrow(location = "bl", which_north = "true",
style = north_arrow_fancy_orienteering)
# Make interactive through plotly? https://plot.ly/ggplot2/maps-sf/
# Using leaflet ----
# Useful help for leaflet
# https://rstudio.github.io/leaflet/choropleths.html
library(htmltools)
# We need to transform the coordinate reference system for help see
# https://mgimond.github.io/Spatial/coordinate-systems-in-r.html
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
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)
# Join with output area data
saveOAprofiles <- left_join(saveOAdata,sf_data, by = c("oaCode" = "oa11cd")) %>%
rename("oa11cd" = "oaCode") %>%
st_as_sf() # reset as sf class
saveLocalAuthList <- unique(saveOAprofiles$ladnm)
# further filter to reduce size of data - no longer used
sf_data2 <- sf_data %>%
dplyr::filter(ladnm %in% saveLocalAuthList)
# create popup first (requires htmltools)
sf_data2$popup_text <-
paste("Output area code: ","<b>", sf_data2$oa11cd, "</b>",
'<br/>', 'Local authority: ', '<b>', sf_data2$ladnm, '</b>', ' ') %>%
lapply(htmltools::HTML)
leaflet(sf_data2) %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addPolygons(color = "blue", fillColor = "blue", fillOpacity = 0.2, weight = 1.5, popup = ~(oa11cd), # popups clicked
label = ~(popup_text), # define labels
labelOptions = labelOptions( # label options
style = list("font-weight" = "normal", padding = "2px 2px"),
direction = "auto"),
highlight = highlightOptions(
weight = 5,
color = "#666",
fillOpacity = 0.7,
bringToFront = TRUE))
dev.off()
# create popup first (requires htmltools)
saveOAprofiles$popup_text <-
paste("Output area code: ","<b>", saveOAprofiles$oa11cd, "</b>",
'<br/>', 'Local authority: ', '<b>', saveOAprofiles$ladnm, '</b>',
'<br/>', 'Peak hours demand (OA): ', '<b>', round(saveOAprofiles$mean4to8), '</b>',' kW',
'<br/>', 'Peak hours demand (Household): ', '<b>', round(saveOAprofiles$mean4to8HH,3), '</b>',' kW') %>%
lapply(htmltools::HTML)
# pal = colorQuantile("Reds", domain = saveOAprofiles$mean4to8, 5)
pal = colorBin("Reds", domain = saveOAprofiles$mean4to8, 5, pretty = TRUE)
leaflet(saveOAprofiles) %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addPolygons(color = "blue", fillColor = ~pal(mean4to8), fillOpacity = 0.9, weight = 1, popup = ~(oa11cd), # popups clicked
label = ~(popup_text), # define labels
labelOptions = labelOptions( # label options
style = list("font-weight" = "normal", padding = "2px 2px"),
direction = "auto"),
highlight = highlightOptions(
weight = 5,
color = "#666",
fillOpacity = 0.7,
bringToFront = TRUE)) %>%
addLegend(pal = pal, values = ~mean4to8HH) %>%
addMiniMap()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment