diff --git a/howTo/gisExample.R b/howTo/gisExample.R
new file mode 100644
index 0000000000000000000000000000000000000000..fadcbb14559c093dcbcb7bc421d85ca2b478a064
--- /dev/null
+++ b/howTo/gisExample.R
@@ -0,0 +1,160 @@
+#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()
+