diff --git a/howTo/openGeogAPI/app3.R b/howTo/openGeogAPI/app3.R
new file mode 100644
index 0000000000000000000000000000000000000000..3d62114e14dea2ce58516720ba20b9c5bd27a975
--- /dev/null
+++ b/howTo/openGeogAPI/app3.R
@@ -0,0 +1,289 @@
+# Load libraries ----
+
+library(readxl)
+library(ggplot2)
+library(sf)
+library(htmltools)
+library(leaflet)
+library(plotly)
+library(dplyr)
+library(reshape2)
+library(RColorBrewer)
+
+# Construct colour palette ----
+industry_pal <- brewer.pal(n = 8, name = "Greys")[4:8] # industry greys, 5 categories
+domestic_pal <- brewer.pal(n = 4, name = "Blues")[2:4] # domestic blues, 3 categories
+transport_pal <- brewer.pal(n = 6, name = "Oranges")[2:6] # transport oranges, 5 categories
+lulucf_pal <- brewer.pal(n = 9, name = "Greens")[4:9] # lulucf greens, 6 categories
+
+# for details
+detailed_pal <- c(industry_pal, domestic_pal, transport_pal, lulucf_pal)
+
+# for totals/outlines
+totals_pal <- c(brewer.pal(n = 9, name = "Greys")[8],
+ brewer.pal(n = 9, name = "Blues")[8],
+ brewer.pal(n = 9, name = "Oranges")[8],
+ brewer.pal(n = 9, name = "Greens")[7:8])
+
+# List local authority areas to load
+# These used to filter emissions data
+# and construct Open Geog API query (geo_query ... to do)
+las_to_load <- c("Southampton","Portsmouth","Winchester",
+ "Eastleigh","Isle of Wight","Fareham",
+ "Gosport","Test Valley","East Hampshire",
+ "Havant","New Forest","Hart","Basingstoke and Deane")
+
+# Load GHG emissions data ----
+
+url_to_get <- "https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/894787/2005-18-uk-local-regional-co2-emissions.xlsx"
+
+tempf <- tempfile(fileext = ".xlsx")
+download.file(url_to_get, tempf, method = "curl")
+
+dt <- readxl::read_xlsx(tempf, sheet = "Full dataset",skip = 1)
+
+x_min <- min(dt$Year)
+x_max <- max(dt$Year)
+
+# Load LA geography ----
+
+# URL as API query - sometimes we don't want all boundaries
+geo_query <- "https://ons-inspire.esriuk.com/arcgis/rest/services/Administrative_Boundaries/Local_Authority_Districts_December_2018_Boundaries_UK_BGC/MapServer/0/query?where=lad18nm%20IN%20(%27Southampton%27,%27Portsmouth%27,%27Winchester%27,%27Eastleigh%27,%27Isle%20of%20Wight%27,%27Fareham%27,%27Gosport%27,%27Test%20Valley%27,%27East%20Hampshire%27,%27Havant%27,%27New%20Forest%27,%27Hart%27,%27Basingstoke%20and%20Deane%27)&outFields=lad18cd,lad18nm,long,lat&outSR=4326&f=geojson"
+
+message("Loading LA geometry from ONS Open Geography API")
+sf_data <- st_read(geo_query)
+#plot(st_geometry(sf_data))
+
+# Useful lookup spatial reference for CRS
+# https://spatialreference.org/ref/epsg/27700/
+st_coord_sys <- st_crs(sf_data) # check coord system
+st_coord_sys # current coord system EPSG: 4326 (is what leaflet wants - good)
+
+# transform the coord system if required
+if(st_coord_sys$epsg != 4326){
+ sf_data <- st_transform(sf_data, "+proj=longlat +datum=WGS84")
+}
+
+# Process data ----
+# Filter local authorities and correct population
+per_capita_dt <- dt %>%
+ filter(Name %in% las_to_load) %>%
+ rename(Population = `Population ('000s, mid-year estimate)`) %>%
+ mutate(Population = Population*1000)
+
+# Add Hampshire totals (sum of all las in las_to_load) to calculate averages across las
+total_hampshire_dt <- dt %>%
+ filter(Name %in% las_to_load) %>%
+ rename(Population = `Population ('000s, mid-year estimate)`) %>%
+ mutate(Population = Population*1000) %>%
+ group_by(Year) %>%
+ summarise_if(is.numeric, sum, na.rm = TRUE) %>%
+ mutate(Name = " Hampshire",
+ `CTRY18NM/RGN18NM` = "South East",
+ `Second Tier Authority` = "Hampshire",
+ Code = "E00000000") %>%
+ select(`CTRY18NM/RGN18NM`,`Second Tier Authority`,Name,Code,Year,everything())
+
+per_capita_dt <- rbind(per_capita_dt,total_hampshire_dt)
+
+per_cap_fun <- function(x) x*1000/per_capita_dt$Population
+per_capita_totals <- data.frame(per_capita_dt[c(1:5,30)], lapply(per_capita_dt[c(11,15,21,28,29,32,33)], per_cap_fun) )
+per_capita_detail <- data.frame(per_capita_dt[c(1:5,30)], lapply(per_capita_dt[c(6:10,12:14,16:20,22:27,29,32,33)], per_cap_fun) )
+# Note that for hampshire LAs there are no emissions in categories Q or S
+
+# Create plot tables ----
+
+pc_detail_plot <- data.frame(per_capita_detail[c(1:5,7:25)])
+pc_detail_plot <- melt(pc_detail_plot, id.vars = c("CTRY18NM.RGN18NM","Second.Tier.Authority","Name","Code","Year"))
+
+pc_general_plot <- data.frame(per_capita_totals[c(3,5,11:13)])
+
+# Calculate emissions and removals for LULUCF separately
+per_capita_lulucf <- pc_detail_plot %>%
+ filter(grepl("LULUCF", variable)) %>%
+ group_by(Name,Year) %>%
+ summarise(
+ LULUCF.emissions = sum(value[value > 0]),
+ LULUCF.removals = sum(value[value < 0]))
+
+per_capita_totals <- left_join(per_capita_totals,per_capita_lulucf, by = c("Name","Year"))
+
+pc_totals_plot <- data.frame(per_capita_totals[c(1:5,7:9,14,15)])
+pc_totals_plot <- melt(pc_totals_plot, id.vars = c("CTRY18NM.RGN18NM","Second.Tier.Authority","Name","Code","Year"))
+
+
+
+# Functions ----
+
+ghg_subset <- function(dt, auth_area = " Hampshire"){
+ dt <- dt %>%
+ filter(dt$Name %in% auth_area)
+ # add filter for categories with input$
+}
+
+
+
+ghg_subset2 <- function(dt, auth_area = " Hampshire", plot_year = 2018){
+ dt <- dt %>%
+ filter(dt$Name %in% auth_area,
+ dt$Year == plot_year)
+ # add filter for categories with input$
+}
+
+# to expand y axis
+expand_ggplot <- function(plot) {
+ y_min <- floor(min(layer_data(plot)$y, na.rm=TRUE))
+ y_max <- ceiling(sum(layer_data(plot)$y, na.rm=TRUE))
+ coord_flip(ylim=c(y_min, y_max))
+ scale_y_continuous(breaks = y_min:y_max)
+}
+
+# Create map (leaflet) ----
+
+# create popup first (requires htmltools)
+library(htmltools)
+sf_data$popup_text <-
+ paste("Locial authority area code: ","<b>", sf_data$lad18cd, "</b>",
+ '<br/>', 'Local authority: ', '<b>', sf_data$lad18nm, '</b>', ' ') %>%
+ lapply(htmltools::HTML)
+
+# App layout ----
+
+ui <- fluidPage(
+
+
+ titlePanel("Hampshire local authority ghg emission visualisation"),
+
+ fluidRow(
+
+ column(5,
+ wellPanel("Click map to select local authority"),
+ leafletOutput("authmap",height = "600px")),
+
+ column(7,
+ wellPanel("Use slider to select year to view in panel below and mouse-over to view more detail"),
+ sliderInput("yearSlider", label = "Year", min = 2005,
+ max = 2018, value = 2018, step = 1, round = TRUE,
+ ticks = FALSE, width = "90%", sep = ""),
+ plotOutput("plot", height = "400px"),
+ plotOutput("bar", height = "150px"))
+ ),
+
+ fluidRow(
+ column(5),
+ column(7)
+
+ )
+
+)
+
+
+ # *Input() functions
+ # *Output() functions
+
+
+
+server <- function(input, output, session) {
+
+ output$authmap <- renderLeaflet({
+
+ leaflet() %>%
+ addTiles() %>% # Add default OpenStreetMap map tiles
+ addPolygons(data = sf_data, layerId = ~(lad18nm),
+ color = "blue", fillColor = "blue", fillOpacity = 0.2, weight = 1.5,
+ #popup = ~(lad18nm), # 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))
+
+ })
+
+ observeEvent(
+
+ input$authmap_shape_click, {
+
+ click <- input$authmap_shape_click
+
+ print(click$id)
+
+ output$selected_auth <- renderText({
+ print(paste0("Local authority area selected: ",click$id))
+ })
+
+ }
+ )
+
+
+ output$plot <- renderPlot({
+
+ plotCaption = paste0("Emissions data: Department for Business, Energy & Industrial Strategy",
+ "\nLocal authority boundary data (2018): ONS Open Geography Portal",
+ "\nVisualisation: rushby.shinyapps.io/LAemissions")
+
+ if(is.null(input$authmap_shape_click$id)){
+ auth_area <- " Hampshire"
+ } else {
+ auth_area <- input$authmap_shape_click$id
+ }
+
+ plotTitle = paste0("Greenhouse gas emissions by sector for ",auth_area, " - ",x_min, " to ",x_max)
+
+ ggplot() +
+ geom_col(ghg_subset(dt = pc_detail_plot, auth_area = auth_area), mapping = aes(x = Year, y = value, fill = variable), position = "stack") +
+ geom_col(ghg_subset(dt = pc_totals_plot, auth_area = auth_area), mapping = aes(x = Year, y = value, colour = variable), fill = "none", position = "stack") +
+ geom_hline(yintercept=0, lwd=0.4, colour="black", linetype = "dashed") +
+ coord_cartesian(ylim = c(-1,15)) +
+ scale_y_continuous(labels=abs) +
+ scale_x_continuous(breaks = 2005:2018) +
+ scale_color_manual(values = totals_pal, guide = FALSE) +
+ scale_fill_manual(values = detailed_pal) +
+ labs(x = "Year",
+ y = "Emissions per capita, tCO2",
+ fill = "Source",
+ title = plotTitle,
+ caption = plotCaption) +
+ theme(legend.position = "right") +
+ theme_classic()
+
+ })
+
+ output$bar <- renderPlot({
+
+ if(is.null(input$authmap_shape_click$id)){
+ auth_area <- " Hampshire"
+ } else {
+ auth_area <- input$authmap_shape_click$id
+ }
+
+ plotTitle = paste0("Greenhouse gas emissions by sector for ",auth_area, " - ",input$yearSlider)
+
+ plot <- ggplot() +
+ geom_col(data = ghg_subset2(dt = pc_detail_plot, auth_area = auth_area, plot_year = input$yearSlider),
+ mapping = aes(x = Name, y = value, fill = variable), position = "stack") +
+ geom_hline(yintercept=0, lwd=0.4, colour="black", linetype = "dashed") +
+ coord_flip() +
+ scale_color_manual(values = totals_pal) +
+ scale_fill_manual(values = detailed_pal) +
+ theme_minimal() +
+ labs(y = "Emissions per capita (tonnes of carbon dioxide/person)",
+ title = plotTitle) +
+ theme(axis.text.y = element_blank(),
+ axis.title.y = element_blank(),
+ axis.ticks.y = element_blank(),
+ #axis.title.x = element_blank(),
+ panel.grid.major.y = element_blank(),
+ legend.position = "none")
+
+ plot + expand_ggplot(plot)
+
+ })
+}
+
+
+shinyApp(ui = ui, server = server)
\ No newline at end of file