diff --git a/howTo/openGeogAPI/local_auth_ghg_plots.R b/howTo/openGeogAPI/local_auth_ghg_plots.R
index 563284eae9daa7df27283379206ae94ec257e8fe..1fed92585405c38602caa5f1ef361fc7003be6ed 100644
--- a/howTo/openGeogAPI/local_auth_ghg_plots.R
+++ b/howTo/openGeogAPI/local_auth_ghg_plots.R
@@ -2,7 +2,7 @@
 
 
 # Load libraries ----
-
+library(here)
 library(readxl)
 library(ggplot2)
 library(plotly)
@@ -85,12 +85,26 @@ per_capita_dt <- dt %>%
   rename(Population = `Population                                              ('000s, mid-year estimate)`) %>%
   mutate(Population = Population*1000)
 
+total_hampshire_dt <- dt  %>%
+  filter(Name %in% las_to_load) %>%
+  group_by(Year) %>%
+  summarise_if(is.numeric, sum, na.rm = TRUE) %>%
+  mutate(Name = "Hampshire")
+
+
 
 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)) %>%
@@ -101,15 +115,9 @@ per_capita_lulucf <- pc_detail_plot %>%
 
 per_capita_totals <- left_join(per_capita_totals,per_capita_lulucf, by = c("Name","Year"))
 
-# Create plot tables ----
-
 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"))
 
-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)])
 
 ## Subset data - Southampton by default
 auth_area  <- "New Forest"
@@ -121,6 +129,12 @@ ghg_subset <- function(dt, auth_area = "Southampton"){
 }
 
 # Construct plots ----
+i <- 1
+plotNames <- list()
+
+for (auth_area in las_to_load) {
+  
+plotName <- paste0("laemissions_by_year_",gsub("[[:space:]]", "_", auth_area))
 
 plotCaption = paste0("Emissions data: Department for Business, Energy & Industrial Strategy",
                      "\nUK local authority and regional carbon dioxide emissions national statistics: 2005 to 2018",
@@ -148,35 +162,51 @@ ggplot() +
   theme(legend.position = "right") +
   theme_classic()
 
-plotly::ggplotly(plot)
+plotNames[i] <-  plotName
+ggsave(paste0(here::here(),"/howTo/openGeogAPI/plots/",plotName,".png"), dpi = 150, width = 12, height = 6, units = "in")
 
+i <- i+1
 
-# By local authority
+}
 
-plot_year <- 2005
+print(plotNames)
 
-plot_data1 <-  pc_detail_plot %>% filter(Year == plot_year)
-plot_data2 <- pc_totals_plot %>% filter(Year == plot_year)
+plotly::ggplotly(plot)
 
-plotTitle = paste0("Greenhouse gas emissions by sector for Hampshire local authorities: ",plot_year)
 
-ggplot() +
-  geom_col(data = plot_data1, aes(x = Name, y = value, fill = variable), position = "stack") +
-  geom_col(data = plot_data2, aes(x = Name, y = value, colour = variable), fill = "none", position = "stack") +
-  geom_hline(yintercept=0, lwd=0.4, colour="black", linetype = "dashed") +
-  coord_cartesian(xlim = c(-1,8)) +
-  #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) +
-  coord_flip() +
-  labs(x = "Local authority",
-       y = "Emissions per capita, tCO2",
-       fill = "Source",
-       title = plotTitle,
-       caption = plotCaption) +
-  theme(legend.position = "none") +
-  theme_classic()
+# By local authority
+years_to_plot <- c(2005,2018)
+
+for(plot_year in years_to_plot) {
+  
+  plotName <- paste0("laemissions_by_auth_",gsub("[[:space:]]", "_", plot_year))
+  
+  
+  plot_data1 <-  pc_detail_plot %>% filter(Year == plot_year)
+  plot_data2 <- pc_totals_plot %>% filter(Year == plot_year)
+  
+  plotTitle = paste0("Greenhouse gas emissions by sector for Hampshire local authorities: ",plot_year)
+  
+  ggplot() +
+    geom_col(data = plot_data1, aes(x = reorder(Name, desc(Name)), y = value, fill = variable), position = "stack") +
+    geom_col(data = plot_data2, aes(x = reorder(Name, desc(Name)), y = value, colour = variable), fill = "none", position = "stack") +
+    geom_hline(yintercept=0, lwd=0.4, colour="black", linetype = "dashed") +
+    coord_flip(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 = "Local authority",
+         y = "Emissions per capita, tCO2",
+         fill = "Source",
+         title = plotTitle,
+         caption = plotCaption) +
+    theme(legend.position = "none") +
+    theme_classic()
+  
+  ggsave(paste0(here::here(),"/howTo/openGeogAPI/plots/",plotName,".png"), dpi = 150, width = 12, height = 6, units = "in")
+  
+}
 
 ggplot() +
   geom_col(data = plot_data1, aes(x = Name, y = value, fill = variable), position = "stack") +